Victor Court
email: victor.court@ifpen.fr

Fernanda Guedes¹

Alexandre Szklo¹

Pedro Rochedo¹

Frédéric Lantz²

Leticia Magalar¹

Eveline Maria

Vásquez Arroyo¹

Climate-energy-water nexus in Brazilian oil refineries

Keywords: Climate-energy-water nexus, oil refineries, Brazil

La collection "Les Cahiers de l'Économie" a pour objectif de présenter les travaux réalisés à IFP Energies nouvelles et IFP School qui traitent d'économie, de finance ou de gestion de la transition énergétique. La forme et le fond peuvent encore être provisoires, notamment pour susciter des échanges de points de vue sur les sujets abordés. Les opinions exprimées dans cette collection appartiennent à leurs auteurs et ne reflètent pas nécessairement le point de vue d'IFP Energies nouvelles ou d'IFP School. Ni ces institutions ni les auteurs n'acceptent une quelconque responsabilité pour les pertes ou dommages éventuellement subis suite à l'utilisation ou à la confiance accordée au contenu de ces publications.

Pour toute information sur le contenu, contacter directement l'auteur.

Introduction

Two of the UN Sustainable Development Goals (SDGs) focus on achieving physical and economic access to energy and water in quantity and quality. SDG7 aims to provide affordable, secure, sustainable and modern energy for all; furthermore, SDG6 aims to provide available and sustainable management of water and sanitation for all (UN, 2016). Energy and water are key elements closely linked to all other sectors within an economy. They also interact closely in many aspects (BIGGS et al., 2015). In this way, achieving the goal of a natural resource will influence the fulfillment of the other goal. For instance, water is needed at all stages of energy production, while water management, treatment and transportation require energy. Moreover, global climate change can add a significant amount of uncertainty to these complex interrelations. Changes in climate variables, such as precipitation and temperature, can affect water and energy resources, increasing their vulnerabilities. Also, the strategies to tackle climate change by reducing (mitigating) greenhouse-gas (GHG) emissions can affect the water-energy nexus (HOWELLS et al., 2013).

For instance, coal-fired thermoelectric plants need water resources, mainly for cooling processes. For these plants, a promising GHG mitigation option could be the installation of amine-based carbon capture (CC) systems [START_REF] Rochedo | Economic analysis under uncertainty of coal fired captureready power plants[END_REF]. However, CC would also increase both water withdrawal and consumption by the thermoelectric plant by more than 100%, which may intensify its vulnerability and affect the water supply to other users downstream from the power plant [START_REF] Zhai | Carbon capture effects on water use at pulverized coal power plants[END_REF]MERSCHMANN et al., 2012). In the case of the production of liquid biofuels, the nexus goes beyond the energy conversion facility, which may also be affected by CO 2 mitigation options, and mostly refers to the biomass production, which usually represents a significant share of water consumption (irrigation) in countries such as Brazil (IEA, 2016). Interestingly enough, the increase of biomass productivity arising from irrigation is an emblematic case of the tradeoff between GHG mitigation and water.

At the end, given all these complex and interconnected relationships, an integrated analysis is needed to evaluate the nexus between energy-water under the challenges associated with climate change (HOWELLS et al., 2013). In addition, each energy sector needs a proper analysis to quantify this nexus. For this study, this analysis is performed at both country and local level. On one hand, the country level provides the basic answer for the primary research question of this study, which is: do CO 2 mitigation options affect the water consumption of an oil refinery system (or even: what could the nexus be between the carbon mitigation cost curve and the water consumption in refineries)? On the other hand, the detailed local level analysis, whose focus is on a specific oil refinery, allows the answering of the secondary question of this study, which is: do the impacts of climate mitigation options on water consumption affect the water supply-demand balance of an oil refinery? Only local level analyses can solve this secondary question, since it requires the proper evaluation of water sources (water supply) and sinks (water users).

In fact, oil refining is an energy-intensive activity, whose greenhouse gas (GHG) emissions are closely related to the combustion and chemical conversion of fossil fuels. The fuel combustion in oil refineries is related to the generation of direct heat, process steam and even electricity, all in stationary sources. The refineries´ technological schemes are complex [START_REF] Gomes | The impact of CO2 taxation on the configuration of new refineries: An application to Brazil[END_REF][START_REF] Coelho | Dealing with petroleum surpluses in Brazil through optimization refining model[END_REF], depending on the characteristics of the feedstocks, the units' capacities, the production profile of the oil products (quantities and specifications), and the choice of technologies to be used [START_REF] Castelo | Abatement costs of CO2 emissions in the Brazilian oil refining sector[END_REF]. For instance, refineries that process heavy crude oils to output light products present process schemes that use more final energy, and in turn emit more GHG (EPA, 2010). In addition, the more stringent the oil derivative specifications, the greater the energy and water consumption of the refining process, due to the need of severe hydro-treatment units, which use the hydrogen produced in units emitting CO 2 from the steam reforming of light hydrocarbons [START_REF] Szklo | Fuel specification, energy consumption and CO2 emission in oil refineries[END_REF][START_REF] Concawe | Refinery energy systems and efficiency[END_REF][START_REF] Sun | Estimation of U.S. refinery water consumption and allocation to refinery products[END_REF]. In 2012, oil refineries accounted for 2.7% of US, 3.2% of European Union and 2.0% of Brazil CO 2 emissions (MCTI, 2013;PETROBRAS, 2013;EPA, 2014).

Nevertheless, the vast majority of the research associated with the nexus between energy and climate in oil refineries has focused on the trade-off between fuel specifications and CO 2 emissions [START_REF] Concawe | Carbon, Sulphur and hydrogen in oil refineries[END_REF][START_REF] Chan | Quantificação e redução de emissões de gases de efeito estufa em uma refinaria de petróleo[END_REF][START_REF] Szklo | Fuel specification, energy consumption and CO2 emission in oil refineries[END_REF][START_REF] Johansson | Assessment of strategies for CO2 abatement in the European petroleum refining industry[END_REF][START_REF] Concawe | Refinery energy systems and efficiency[END_REF]. In the case of the nexus between energy and water, there are some studies on the relationship between water consumption and energy use in oil refineries [START_REF] Hightower | The energy challenge[END_REF]IPIECA, 2010;[START_REF] Hwang | Water network synthesis in refinery[END_REF][START_REF] Pan | A supply chain based assessment of water issues in the coal industry in China[END_REF][START_REF] Mughees | Application of Water Pinch Technology in Minimization of Water Consumption at a Refinery[END_REF][START_REF] Sun | Estimation of U.S. refinery water consumption and allocation to refinery products[END_REF]. Previous research has also focused on the implementation of CO 2 capture in oil refineries and its abatement cost, as [START_REF] Rochedo | Carbon capture potential and costs in Brazil[END_REF] have done, but it has failed to explore the water nexus with CO 2 mitigation options in oil refineries.

At the end, few attempts have been made to quantify the relationship between climate (CO 2 emission mitigation) and water-energy use in oil refineries. Table 1 provides a brief summary of the opportunities for CO 2 abatement measures in oil refineries' processing units, and their likely impact on water consumption. It highlights the signs of the impacts that are quantified later in this study (positive and negative signs) through the use of simulation tools for all Brazilian refineries and for a specific refinery in detail.

- - - - - VDU - - - - - CRU - - - HDT - - - - - - RFCC - - FCC - - + DCU - - -

HGU +

Note: The positive sign means that the measure has a positive impact on water consumption, ie it helps to reduce this consumption.

On the other hand, the negative sign has a negative impact, indicating an increase in water consumption with the application of the measure. Regarding water use, oil refining requires considerable amounts of water, which vary significantly between refineries, depending on the process configuration [START_REF] Wu | Consumptive water use in the production of ethanol and petroleum gasoline[END_REF][START_REF] Sun | Estimation of U.S. refinery water consumption and allocation to refinery products[END_REF], the petroleum specification (e.g., density, sulphur content, total acid number) [START_REF] Sun | Estimation of U.S. refinery water consumption and allocation to refinery products[END_REF], and the products' requirements. Within this context, in Brazilian refineries, water requirements deserve attention due to the processing of heavy-to-medium crude oils, as well as to the increasingly stringent specifications of fuels that require the implementation of hydrotreatment units, associated with the water-intensive steam reforming process [START_REF] Castelo | CO2 emissions abatement costs of reducing natural gas flaring in Brazil by investing in offshore GTL plants producing premium diesel[END_REF][START_REF] Szklo | Fundamentos do Refino do Petróleo[END_REF][START_REF] Barros | Petroleum refining flexibility and cost to address the risk of ethanol supply disruptions: The case of Brazil[END_REF] 1 .

ADU

In addition, in Brazil, water availability is not evenly distributed. While the northern region holds more than 80% of all water availability, the basins located in large urban centers, in the with the objective of meeting the demand for gasoline and fuel oil in major urban centers (also close to Brazil's coast). However, due to the increasing diesel demand after the 1980s (BORBA et al., 2017), as well as the ramp up of medium-to-heavy crude oil production in Brazilian offshore basins in the 1980s and 1990s [START_REF] Hallack | Curve-fitting variants to model Brazil's crude oil offshore post-salt production[END_REF], the refining schemes of existing refineries were altered to convert the heaviest fractions of crude into medium cutse.g., by adding delayed coking units and severe hydrotreatment processes (which remove nitrogen compounds, high sulfur compounds and aromatic rings) [START_REF] Szklo | Fuel specification, energy consumption and CO2 emission in oil refineries[END_REF][START_REF] Szklo | Fundamentos do Refino do Petróleo[END_REF].

for 54% of Brazil's refining capacity, are located in areas classified as having critical water availability. To sum up, oil refineries are major CO 2 emitters and are usually located in sites with already critical availability supply of water. This is the case for Brazil, but it can be seen worldwide [START_REF] Oil | World Wide Refining Survey[END_REF]. Some CO 2 emission mitigation measures can positively or negatively influence both withdrawal and water consumption at refineries. Within this context, this study aims to quantify the extent to which these measures can impact refineries' water use.

First, by developing an energy, CO 2 , H 2 , water balance simulator for Brazilian oil refineries, and applying it to different scenarios of CO 2 mitigation, this study evaluates the CO 2 mitigationenergy-water nexus at the country level. Then, this study analyzes the water supply of the hydrographic basin in which the largest Brazilian refinery, REPLAN, is located and how it behaves over time. This detailed analysis, at a local level, not only quantifies the water withdrawal impacts of CO 2 mitigation options, but also identifies whether or not these impacts could be overcome by the current water supply of this specific refinery.

The next section presents the methodology used to carry out the analysis, as well as a brief description of the simulation tools applied. Section 3 discusses the results obtained. Lastly, the final remarks of the study, highlighting also its limitations, are presented.

Methods

Methodological Procedure

The methodological procedure applied by this study consisted of the following steps. The first step is the definition of GHG mitigation measures in oil refineries, according to the scientific literature and the experience of the authors regarding Brazil's oil refineries 2 . Then, the study simulates the baseline case for estimating the scenario without CO 2 mitigation, and simulates the impacts of introducing CO 2 price scenarios into oil refineries, in terms of CO 2 emissions, final energy use and water consumption. This provides the carbon mitigation cost curve, and also helps to identify (quantify) the impacts of GHG mitigation options on the water consumption for all Brazilian refineries. However, this step is not able to detail the water supply balance at a local level. Therefore, using the mitigation cost curve for Brazilian refineries (mentioned above), the proposed procedure includes a last step for detailing the case of the largest Brazilian oil refinery, REPLAN. This step not only quantifies the water consumption impacts of CO 2 mitigation options, but also identifies whether these impacts could be overcome by the current water supply of REPLAN. In summary, the steps include:

1. To assess the CO 2 mitigation options focusing on the saving potential for the consumption of fuels, steam, electricity and H 2 .

2. To estimate the energy and mass balances for a baseline case, including final energy consumption, CO 2 emissions and water requirements. This case does not consider the application of CO 2 emission mitigation options (e.g., fuel switch, fuel saving and carbon capture). This study applies an energy and mass balance simulatorthe so-called 2 See, for instance, GUEDES (2015).

"CAESAR -Carbon and Energy Strategy for Refineries" tool, which is briefly described here, and better described in the Supplementary Material.

3. To run CAESAR with CO 2 emission prices 3 of 25, 50, 100 and 200 US$/tCO 2 . In this case, the CO 2 emissions mitigation options are selected according to their marginal abatement costthat is, technological options with costs lower than or equal to the exogenously established CO 2 price are automatically selected by the simulation tool, allowing the construction of a CO 2 average abatement cost curve for all Brazilian oil refineries. The Supplementary Material provides the basic equation associated with the estimation of the abatement cost.

4. To develop a case study for REPLAN, the largest Brazilian refinery in terms of processing capacity, thus quantifying the water stress in detail, or locally. This allows investigating whether mitigation measures that were selected in step 3 can be adopted in cases where a greater water withdrawal is required. This case study is performed using the software tool Water Evaluation and Planning -WEAP (see section 2.2. and Supplementary Material).

CAESAR tool -Carbon and Energy Strategy Analysis for

Refineries

The tool used for evaluating all Brazilian refineries, without detailing the water supply-demand balance at a local level, is the simulator CAESAR -Carbon and Energy Strategy Analysis for Refineries. It was originally developed by [START_REF] Tolmasquim | Matriz Brasileira na Virada do Milênio[END_REF], later being used by the Brazilian Government in its Long-Term Energy Plan 2030 (EPE, 2007). Finally, it was updated by [START_REF] Guedes | Avaliação de alternativas para redução do uso final de energia no setor de refino de petróleo brasileiro e estimativa de custos de abatimento de emissões de gases de efeito estufa[END_REF] and by [START_REF] Vásquez | Incorporação do nexo energia-água em um modelo de otimização da expansão do sistema energético brasileiro[END_REF] and [START_REF] Magalar | Análise da disponibilidade hídrica como um fator de restrição à produção de derivados em refinarias[END_REF] for incorporating water balances.

The simulation is performed within Excel (visual basic), and relies on refining schemes, including the following units' energy and mass balances: atmospheric distillation, vacuum distillation, alkylation, atmospheric residue delayed coking, vacuum residue delayed coking, propane desasphalter, catalytic reformer, fluid catalytic cracker, hydrocracker, residue fluid catalytic cracker, hydrotreaters (naphtha, diesel, kerosene and instable products), hydrotreatment of finished gasoline, lube unit, and hydrogen generation unit. The processing units' capacities are determined, as well as the processed feedstocks, specific utilities consumption (steam, fuel and hydrogen) and specific water consumption. The outputs of the tool consist of the final energy consumption, CO 2 emissions, oil product output, and refineries' water consumption and withdrawal.

3 They represent an established price to be paid for a given amount of CO 2 emitted.

Therefore, CAESAR is a bottom-up model mostly based on the simulation of the mass (water, H 2 ) and energy balances of Brazilian oil refineries. It has an additional feature for optimizing the energy consumption aimed at minimizing the cost of operation of oil refineries. The model also includes a list of CO 2 mitigation options, which are detailed according to the processing units in which they can be implemented, their potential for saving fuel and/or electricity, their investment, operation and maintenance costs, and their penetration rates. In total, 204 options of technologies are available in CAESAR (see Supplementary Material for detailed data).

For the carbon price scenarios, CO 2 emission prices were exogenously introduced into the simulator, which also affected the optimization problem that finds the least-cost fuel mix of

WEAP

Before using the tool WEAP, REPLAN mass and energy balances were simulated in the abovedescribed tool, CAESAR. This aimed to quantify the impacts of CO 2 emission mitigation options on the water required by REPLAN. Then, the results of water withdrawals obtained in CAESAR were inserted as input into the WEAP tool. This is a tool for integrated water resources management (IWRM) developed by the Stockholm Environmental Institute (SEI).

WEAP integrates physical hydrological processes with water withdrawal management and infrastructure, as well as environmental and economic aspects of water planning. Its simulations are based on scenarios that can be analyzed according to different trends in hydrology, water use and demand, demography, technology, operating rules and water management policies [START_REF] Sieber | Recent advances in the science and technology of ultra-low sulfur diesel (ULSD) production[END_REF].

The WEAP analysis consists of, firstly, configuring the time horizon, catchment areas, system components and configuration of the problem to be evaluated. Then, the model is used to simulate alternative scenarios to assess the impact of different water supply and demand management options, as well as evaluate the water availability within a region of study.

The model simulates the use of water in hydrological basins by using a linear programming algorithm, which aims to maximize the water delivered to demand sites, according to a set of priorities defined by the user. When water is limited, the algorithm is formulated to progressively constrain water allocation to the lowest priority demand sites. More details of the model can be found in SIEBER and PURKEY, 2015. See the Supplementary Material for further details on how WEAP is calibrated and used by this study.

Input Data

Brazilian Case Study in CAESAR

The analysis performed by this study was based on the current Brazilian oil refinery system, thus, no greenfield refinery was constructed in the simulation. The mass and energy balances rely on the breakdown in processing units, which have specific characteristics. The capacity of these units is shown in Table 2. The average utilization factor of the atmospheric distillation unit was set as 70%, following MME (2018). Source: OIL AND GAS JOURNAL (2018)

Table 3 shows the estimates for Brazilian refineries' typical utility consumption (negative values mean a net production of the utility by the unit). Although there are variations in the specific energy consumption of utilities for the same unit, depending on the supplier of the technology, local characteristics or even different design considerations, the values adopted in CAESAR seek to represent a typical Brazilian unit. For Brazilian oil refineries, as of 2017, coefficients of water withdrawals per process unit were determined, as indicated in Table 4. The coefficients consist of low-pressure steam (LP Steam), medium-pressure steam (MP Steam) and high-pressure steam (HP Steam), related to the process units. The water balance also includes the water consumed in the cooling system (CWcooling water) and the volume of demineralized water used in the boiler (BFWboiler feed water) per barrel of oil processed. Figure 2 shows the basic water balance applied in the simulation tool. After obtaining the steam demand, CW and BFW, parameters were adopted to estimate the water consumed by the refinery, based on ANZE (2013). They are composed of the make-up water used in the cooling towers, the make-up water for the boilers and the water used by the processes (water incorporated into products, for instance, in the production of H 2 , because of steam reforming and water gas shift, during chemical reactions). A value of 1.7% was considered for the cooling system's total circulating water, according to typical Brazilian oil refineries' concentration ratios [START_REF] Magalar | Análise da disponibilidade hídrica como um fator de restrição à produção de derivados em refinarias[END_REF]. For the boiler water make-up, a value of 49.7% was applied to the sum of the amount of water used in the boilers (BFW) and the total amount of steam consumed in process units. Steam consumed is defined as lost steam that did not return as condensate. For this, a value of 33% of all generated steam was used (VÁSQUEZ [START_REF] Vásquez | A DISPONIBILIDADE HÍDRICA COMO POSSÍVEL RESTRIÇÃO AMBIENTAL PARA A PRODUÇÃO DE DERIVADOS DE PETRÓLEO NO BRASIL[END_REF]. Equation ( 1) summarizes these assumptions and the water balance.

( ∑ ) ∑ [ ∑ ] (1)
Where "CW i " represents the cooling system's circulating water of each process unit "i";

"BFW i " is the amount of water used in boilers of each process unit "i"; "Steam i " is steam consumed in each process unit "i".

REPLAN Case Study in WEAP

REPLAN is the largest Brazilian refinery in terms of processing capacity (66 thousand m³/day) (ANP, 2018). This refinery is located in Paulínia in the state of São Paulo and is placed in the Piracicaba, Capivari and Jundiaí River Basin (PCJ), which is classified as "critical" in relation to water availability [START_REF] Magalar | Análise da disponibilidade hídrica como um fator de restrição à produção de derivados em refinarias[END_REF].

In this study, the water availability of Jaguari basin was calculated by simulating a water balance between the inflows and outflows of its drainage area over time. The Jaguari River basin was chosen due to the catchment point of REPLAN being located in this river. In addition, the basin of the Camanduacaia, Ribeirão do Pinhal rivers was integrated into the case study because these rivers are tributaries of the Jaguari River.

The water balance method chosen by this study was the simplified coefficient Method -Rainfall Runoff, in which water requirements are calculated based on evapotranspiration and precipitation data. Twenty-six rainfall stations were evaluated within the three catchment areas:

Jaguari River catchment and its tributary rivers, Camanducaia and Pinhal. The data loaded after treatment of the missing data and outliers was the monthly average rainfall. In order to use the mean value of evapotranspiration for each catchment area, the monthly average of all municipalities in each area was calculated. For more details of data, see Supplementary

Material.

The water outputs considered in this study are the projected demands for public supply, industry, irrigation and for animal husbandry. These demands were identified and projected to the year 2040 to assess the extent to which water availability changes as a function of the multiple uses of water within the PCJ basin and whether the REPLAN could be impacted.

The demand for water for urban supply was calculated using a coefficient of water demand per inhabitant per day that was adjusted to account for the water losses in distribution. The same coefficient was used for the projection of water demand for future public supply. The method used for the estimation of the population of each city is described in the Supplementary Material.

Water consumption for animal husbandry was calculated from data on the number of animals per city and then calculated the product of the effective number of herds by a per capita coefficient of daily water consumption known as equivalent cattle for water demand. In order to estimate the industrial demand, the volume of water granted by industry in the water agency was consulted.

The demand for irrigation is calculated by multiplying the area under cultivation by the difference between the water requirement of the crop and the precipitation occurring over the cultivated area. For this, it is necessary to know the water demand of each crop, which is calculated from the reference evapotranspiration and crop coefficient.

Results

Baseline Scenario for all Brazilian Refineries

The total consumption of utilities and fuels in existing Brazilian refineries is shown in Table 5.

Negative values indicate exports or utility surpluses, while positive values indicate consumption of utilities. From the utilities consumption, it was possible to determine the fuel consumption. The refinery fuels include natural gas, refinery gas, fuel oil, naphtha and petcoke. Electricity purchased from the grid was also accounted for, either from those refineries that do not have cogeneration or from the excess demand in relation to the capacity of cogeneration units. Natural gas is used for producing hydrogen in HGUs, electricity in cogeneration units, and steam in boilers and direct heating in process units. Refinery gas and fuel oil were accounted for direct heating in process units. In general, leftover refinery gas was directed toward flare emissions accounting.

Furthermore, a 100% flare combustion efficiency was assumed to be conservative on the GHG emission estimates. Finally, the consumption of petcoke was accounted for in FCC and RFCC units. Table 6 shows the estimation of the final energy consumption for the existing Brazilian refineries. By dividing the total emissions by the processed feed, this study estimated an emission intensity of 0.4 tCO 2 /t oil, which is compatible with the 2012 data presented by the Brazilian oil company that owns most of the country's refineries (PETROBRAS, 2013) 6 . Just for comparison, worldwide several works in the literature present CO 2 emission intensities of oil refineries hovering between 0.1 and 0.4 tCO 2 /t of oil processed, with an average of 0.22 [START_REF] Concawe | Impact of product quality and demand evolution on EU refineries at the 2020 horizon CO2 emissions trend and mitigation options[END_REF]IEAGHG, 2008;[START_REF] Straelen | [END_REF][START_REF] Dnv | Global Technology Roadmap for CCS in Industry -Sectoral Assessment: Refineries[END_REF]. For example, the US has an average emission of 0.33 tCO2/t of processed oil, while the European Union has an average value of 0.27 (EPA, 2014).

Finally, concerning the relationship between CO 2 emissions and water withdrawals, the estimative for the baseline scenario is 0.62 tCO 2 /m³. Figure 3 and Figure 4 present, for this scenario, the most representative units in terms of water consumption and CO 2 emissions, respectively.

6 Equal to 0.45 tCO 2 / t of oil processed in 2012. Of course, this intensity may vary slightly among years given the focus of the ADU campaign (in our study we focused on diesel), the possible maintenance of downstream units, which can affect the utilization factor of oil refineries (we used the ADU average utilization factor of 2017, equal to 70%), and the crudes processed in the refineries. In our study, we have considered the ramp-up of a lighter and sweeter feed that has been made available in Brazil in the last five years, from pre-salt fields. That is why we run the model with 4% of the feed from paraffinic oils from Saudi Arabia; 2% from ultra-light African crudes; 32% from Brazilian heavy crudes, and the remaining 62% from medium-to-slightly light Brazilian crudes, mostly from pre-salt fields. Therefore, the feedstock blend has become lighter than it was in 2012. According to [START_REF] Szklo | Fuel specification, energy consumption and CO2 emission in oil refineries[END_REF], most CO 2 emissions from Brazilian refineries come from burning fuels. Interestingly, the fuel consumption of refineries in absolute terms concentrates on few processes, which are not the most energy intensive (in terms of energy consumption per barrel) but process large volumes of feedstock. Typically, atmospheric and vacuum distillation units account for 35-40% of a refinery's final energy use (API, 2000) because any barrel of oil entering a refinery passes through the topping separation units. This explains their share of CO 2 emissions. Also, for global refining, between 16% and 20% of the total are non-energy emissions associated with the chemical reactions of hydrogen production and cracking of the FCC [START_REF] Szklo | Fuel specification, energy consumption and CO2 emission in oil refineries[END_REF]. This average figure agrees with our findings for Brazil. Finally, severe hydrotreatment (for unstable and unfinished distillates) results in both higher water consumption and CO 2 emissions due to the severity (temperature higher than 450°C, H 2 partial pressure up to 21 MPa, and low liquid hourly space velocity7 and hydrogen pressure) under which reactions must happen [START_REF] Gary | Petroleum Refining: Technology and Economics[END_REF]STANISLAU et al, 2010).

CO 2 Price Scenarios for all Brazilian Refineries

As described above, four CO 2 price scenarios were simulated in CAESAR. According to the levelized cost of mitigation options on the database of the tool, different options were selected for each scenario. Moreover, the fuel mix also changed to minimize operational costs considering the CO 2 prices (and the emission factors of each possible fuel to be used). Table 9 summarizes the results for different CO 2 prices, and Figure 5 shows CO 2 emissions and water requirements for different CO 2 emission prices scenarios. The most significant CO 2 emission abatement occurs at 25 and 100 US$/tCO 2 , 10% and 26%, respectively, compared to the baseline. This is explained by the total abatement potential of the technologies found in the cost ranges 0-25 US$/tCO 2 and 50-100 US$/tCO 2 , equal to 143.5 and 205.7 MtCO 2 (see Supplementary Material). In respect to water requirements, a slight change of less than 1% occurs between the baseline scenario and 25 US$/tCO 2 scenario. In other scenarios, the water withdrawals remain practically stable, with a small change, less than 0.5% in the 100US$/tCO 2 scenario. To better illustrate the relationship between the abatement costs and the accumulated abatement potential, the abatement cost curve (Figure 6) was produced, including the 204 technologies considered in the study. In the end, the findings of this study show that the co-benefits of GHG abatement measures that also reduce steam consumption (e.g., reduction of heat storage between ADU and VDU, steam fouling reduction in ADU, installation of vacuum pumps to replace steam injectors in ADU, increase AGR solvent concentration in HDS G, replace steam drive for electric in HDT N, and installation of CO-kiln in regenerative tower HRSG in FCC), which were chosen 8 by our simulations, were offset by the water consumption increase related to CC options, especially in HGU. In summary, at a national level and on average, CO 2 mitigation impacts on water use by oil refineries in Brazil are neutral. Figure 7 illustrates how steam consumption reduction from some mitigation measures is overcome by the increase required with CC implementation.

Case Study: REPLAN

At a local level, for the largest Brazilian oil refinery, the water balance undertaken showed that, although there was no unmet water demand at the REPLAN's catchment point, the conflict between the multiple water users in the basin should intensify. This is due to the trend in the river flow being progressively closer to the critical threshold of 50% of the minimum availability (Q 7,10 ). In addition, it was observed that the point of flow observation at Jaguari River faces instants when the flow must be restricted. This means that REPLAN may sometimes suffer impacts on its operation due to a 30% reduction in the volume of water it receives from the Jaguari River. Figure 8 shows an increase in water withdrawal in all scenarios. Although some mitigation measures reduce steam consumption, in the refinery's overall water balance, this reduction is offset by the increase in the demand for boiler feed water and for cooling. Mitigation measures costing up to US$ 25/MtCO 2 were the ones that most demanded water due to the increase in boiler feed water need, which was 1.4% more than in the baseline scenario. In addition, the slight increase that occurred between scenarios US$ 50/MtCO 2 and 100 was due to the implementation of CC, which increased the demand for cooling water.

Nevertheless, as a final balance, a reduction of less than 1% was obtained when implementing all CO 2 mitigation measures in REPLAN. This means that, contrary to what happens in other energy sectors [START_REF] Zhai | Carbon capture effects on water use at pulverized coal power plants[END_REF]MERSCHMANN et al., 2012), the implementation of CO 2 abatement in oil refineries has no significant impact on water consumption (no negative trade-off). However, this also means that the water stress in oil refineries should be dealt with measures not directly linked to CO 2 abatement (no significant co-benefits). This is valid both at local and country levels.

Final Remarks

This study developed an energy, CO 2 , H 2 , water balance simulator for Brazilian oil refineries, and applied it to different scenarios of CO 2 mitigation (at 25, 50, 100 and 200 US$/tCO 2 )

aiming at investigating the climate-energy-water nexus. A Baseline scenario, i.e., a scenario without CO 2 prices was also elaborated. Results for both scenarios included final energy consumption, CO 2 emissions and water requirements. The most significant reductions in CO 2 emissions were due to the implementation of the carbon capture. However, this option offsets the co-benefits of CO 2 abatement measures that reduced the water requirements of Brazilian oil refineries, especially those already located in areas under water supply stress, such as the largest refinery in Brazil (REPLAN), whose water balance with carbon mitigation options was detailed in this study.

Nevertheless, as this study focused on the impacts of CO 2 mitigation options on water requirements, it was not able to follow the reverse path of the nexus: from climate to water availability. This means that climate change can affect the water availability to oil refineries (water supply, instead of water demand side). Hence, future studies could focus on this issue, also including the analysis of alternatives to regularize river flows to deal with climate impacts on water supply. Another idea could be optimizing refineries for minimizing water consumption (or withdrawals).

It is also worth noting that this study tried to validate the findings of the tools used by comparing them to real data from Brazil. However, an important issue for the simulation tool is to calibrate the feedstock blend to be run, and the focus of the refinery operation. As of today, although the Brazilian refinery system, on average, focuses on diesel optimization (e.g. when establishing the distillation cuts), single refineries can present a different feature (e.g. focusing on lube oils or petrochemicals). Similarly, the yearly focus of the average refinery operation on diesel does not mean that this is valid for all days of the year.

Finally, although the 204 CO 2 mitigation options considered by this study represent an extensive list of measures, there are always new possibilities to be assessed. For example, some studies have evaluated the use of renewable energy sources to supply the energy demand (PINSKE et al., 2012) and the hydrogen consumption (SILVA, 2017) of oil refineries.

  atmospheric distillation unit; VDUvacuum distillation unit; CRUcatalytic reforming unit; HDThydrotreatment unit; RFCCresid fluid catalytic cracking; FCCfluid catalytic cracking; DCUdelayed coking unit; HGUhydrogen generation unit; CC: CO2 capture Source: WORREL and GALITSKY (2003); WORREL and GALITSKY (2005); CONCAWE (2008); BERGH, 2012; MORROW III et al. (2013)

Figure 1 -

 1 Figure 1 -Refineries locations and water criticality indicator of sub-basins in BrazilSource: Based on ANA(2017)

Figure 2 -

 2 Figure 2 -Water balance in CAESAR Source: Based on IPIECA (2010)

  After all climatic parameters, data on land use and water demands are inserted into the model, the observed values of the fluviometric stations are compared with the flow data modeled by WEAP. From the observed and simulated flow data, two calibration indices are calculated, the Nash-Sutcliffe efficiency index and the BIAS index.To evaluate the water availability of the REPLAN catchment area, a minimum ecological flow was defined. The minimum flows most commonly used in Brazil are Q 7where the drainage area is located. According to the water resources committee of the PCJ (CBH-PCJ, 2000), areas considered critical are those in which the total water demand exceeds 50% of the minimum availability Q 7,10. In addition, the water resources policy in the state of Sao Paulo determines that the volume of water withdrawal in the industrial sector should be reduced if the flow of the Jaguari River reaches the minimum flow established in specific gauge stations. Therefore, in the water balance simulation done by this study, the analysis tried to find out if periods of restriction of water for REPLAN could happen.

Figure 3 -

 3 Figure 3 -Water consumption per processing unit in the baseline scenario ADUatmospheric distillation unit; VDUvacuum distillation unit; ; FCCfluid catalytic cracking; RFCCresid fluid catalytic cracking; ALKalkylation unit; DCUdelayed coking unit; CRUcatalytic reforming unit; HDS Ggasoline hydrodesulphurization unit; HDT Nnaphtha hydrotreatment unit; HDT Qkerosene hydrotreatment unit; HDT Ddiesel hydrotreatment unit; HDT Isevere hydrotreatment unit; LUBlubricants unit; HGUhydrogen generation unit

Figure 5 -

 5 Figure 5 -Water withdrawal versus CO 2 emissions

Figure 6 -

 6 Figure 6 -Abatement cost curve for Brazilian refineries Note: HGU -Hydrogen generation unit; FCC -Fluidized catalytic cracking.

Figure 7 -

 7 Figure 7 -Steam requirements impacts of CO 2 mitigation options

Figure 8 -

 8 Figure 8 -Water withdrawal versus CO 2 abatement in REPLAN

Table 1 -Qualitative Impacts on Water Consumption of CO 2 Mitigation Options in Oil Refineries Process Unit Heat Integration

 1

	Reduce Boiler	Blowdown/Water	Treatment	Improved	Maintenance/Steam	Lines & Traps	Reduce Stand-by	Boiler Requirements	Increase Steam Line	Insulation	Recover Blowdown	Steam	CC
	ADU

  refineries. Prices of 25, 50, 100 and 200 US$/tCO 2 were considered, thus building five different scenarios for the current configuration of Brazilian oil refineries. As 204 CO 2 emission

mitigation options are available in the simulator, their abatement costs range from negative values, which represent "non-regret" measures, to values above 100 US$/tCO 2 . The highest cost measures would hardly come into effect without economic incentives or more robust technological learning.

Therefore, depending on the CO 2 emission price applied, the tool automatically selects different GHG mitigation options from the set list available, affecting the final energy use, CO 2 emissions and water consumption. The Supplementary Material includes the basic data of the model and a description of how to run it.

Table 2 -Brazilian Process Unit Capacities as of December 2017 Unit Capacity (barrels/d)

 2

	ADU	2,138,000
	VDU	804,740
	FCC	378,729
	RFCC	123,158
	ALK	6,290
	DCU	115,319
	CRU	2,386
	HDS G	3,054
	HDT N	10,528
	HDT Q	28,125
	HDT D	200,041
	HDT I	11,698
	LUB	20,009
	HGU	126
	ADU -atmospheric distillation unit; VDU -vacuum distillation unit; ; FCC -fluid catalytic cracking; RFCC -resid
	fluid catalytic cracking; ALK -alkylation unit; DCU -delayed coking unit; CRU -catalytic reforming unit; HDS G-gasoline

hydrodesulphurization unit; HDT Nnaphtha hydrotreatment unit; HDT Qkerosene hydrotreatment unit; HDT Ddiesel hydrotreatment unit; HDT Isevere hydrotreatment unit; LUBlubricants unit; HGUhydrogen generation unit (in this case, the capacity is given in MMcfd)

Table 3 -Process Units' Utilities Specific Energy Consumption

 3

	Unit	HP Steam	MP Steam	LP Steam	Electricity	Fuel	Coke	H 2 Consumption	H 2 Production	BFW	CW
		kg/bbl kg/bbl kg/bbl	kWh/bbl	MJ/bbl MJ/bbl	m³/bbl	m³/bbl	m³/bbl m³/bbl
	ADU	0.00	11.00	0.00	0.60	127.00	0.00	0.00	0.00	0.02	0.35
	VDU	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.05	0.35
	FCC -16.00 20.00	-3.60	8.80	0.00	368.00	0.00	0.00	0.07	1.00
	RFC C	-18.00	0.00	0.00	1.00	0.00	368.00	0.00	0.00	0.07	1.00
	ALQ	0.00	90.00	0.00	9.00	0.00	0.00	0.00	0.00	0.05	7.00
	CRU -15.60	0.00	0.00	10.00	382.00	0.00	-48.00	0.00	0.02	1.74
	DCU	0.00	-18.40	0.00	3.60	126.00	0.00	0.00	0.00	0.06	2.03
	HDS G	3.00	0.00	0.00	2.00	105.00	0.00	4.00	0.00	0.04	0.96
	HDT N	3.00	0.00	0.00	2.00	105.00	0.00	7.00	0.00	0.01	0.19
	HDT Q	4.00	0.00	0.00	3.00	158.00	0.00	7.00	0.00	0.18	0.49
	HDT D	4.00	0.00	0.00	3.00	158.00	0.00	7.00	0.00	0.04	0.73
	HDT I	5.00	0.00	0.00	6.00	211.00	0.00	17.00	0.00	0.05	0.71
	LUB	0.00	1.60	5.60	1.60	135.00	0.00	0.00	0.00	0.05	1.00
	HGU	0.00	0.00	0.00	0.00	2.55	0.00	0.00	0.16	0.00	0.00
			HP -high pressure; MP -medium pressure; LP -low pressure; BFW -boiler feed water; CW -cooling
						water
			Source: Based on HYDROCARBON PROCESSING (2008); MEYERS (2004); GARY AND
					HANDWERK (2001); STANISLAUS et al. (2010)

Table 4 -Water Use Coefficients per Process Unit

 4

		CW	BFW	LP Steam	MP Steam	HP Steam
		(m³/bbl)	(m³/bbl)	(kg/bbl)	(kg/bbl)	(kg/bbl)
	ADU	0.3	0.02	-	11.0	-
	VDU	0.3	0.05	-	-	-
	FCC	1.0	0.07	3.6	20.0	16.0
	RFCC	1.0	0.07	-	-	18.0
	ALQ	7.0	0.05	-	90.0	-
	CRU	1.7	0.02	-	-	15.6
	DCU	2.0	0.06	-	18.4	-
	HDS G	1.0	0.04	-	-	3.0
	HDT N	0.2	0.01	-	-	3.0
	HDT Q	0.5	0.18	-	-	4.0
	HDT D	0.7	0.04	-	-	4.0
	HDT I	0.7	0.05	-	-	5.0
	LUB	1.0	0.05	-	-	-
	HGU	-	-	-	-	-

HPhigh pressure; MPmedium pressure; LPlow pressure; BFWboiler feed water; CWcooling water Source: VÁSQUEZ

[START_REF] Vásquez | A DISPONIBILIDADE HÍDRICA COMO POSSÍVEL RESTRIÇÃO AMBIENTAL PARA A PRODUÇÃO DE DERIVADOS DE PETRÓLEO NO BRASIL[END_REF]

Table 5

 5

	-Utilities Consumption
	HP Steam (kt/year)	-752,5
	MP Steam (kt/year)	7,340.4
	LP Steam (kt/year)	-375.9
	Electricity (GWh/year) 12,162.7
	Fuel (TJ/year)	284,018.2
	Coke (TJ/year)	58,952.6
	H 2 (M Nm³/year)	6,116.5
	HP -high pressure; MP -medium pressure; LP -low pressure

Table 6 -Final Energy Consumption -Baseline (PJ/year)

 6

	Natural Gas	367.8

Table 7 -Water Requirements -Baseline

 7 It was also possible to estimate the CO 2 emissions of Brazilian refineries as of 2017, through the multiplication of the emission factors reported byIPCC (2006) of the respective fuels used by Brazilian refineries (Table8). For electricity's CO 2 emissions, the average Brazilian grid emission factor for 2017 was considered, equal to 92.7 tCO 2 /GWh (MCTIC, 2018).

	BFW (t/h)	8,809.2
	CW (t/h)	118,486.2
	Steam (t/h)	3,720.4
	Condensed Steam (t/h)	2,492.7
	BFW spent (t/h)	10,036.9
	BFW Make-up (%)	49.7
	BFW Make-up (t/h)	4,988.4
	CW Make-up (%)	1.7
	CW Make-up (t/h)	2,014.3
	Consumption (t/h)	1,610.5
	Withdrawal (t/h)	7,002.6
	Consumption (m³/bbl)	24.9
	Withdrawal (m³/bbl)	108.2
	Consumption (km³/year)	14107.9
	Withdrawal (km³/year)	61351.7
	BFW -Boiler feed water; CW -Cooling water

Table 8 -CO 2 Emissions (MtCO 2 /year) -Baseline

 8

	Natural Gas	20.6
	Refinery Gas	4.9
	Fuel Oil	6.6
	Coke	5.7
	Grid Electricity	0.4

Table 9 -Summary of Results Final Energy Use (PJ/year) CO 2 Emission Price (US$/tCO 2 )

 9

		Baseline	25	50	100	200
	Natural Gas	367.84	367.84	367.84	367.84	367.84
	Refinery Gas	84.39	84.39	84.39	84.39	84.39
	Fuel oil	85.26	78.76	78.22	71.32	62.05
	Coke	58.95	58.95	58.95	58.95	58.95
	TOTAL	596.40	589.94	589.40	582.50	573.23
	Grid Eletricity (GWh/year)	7252.71 7393.02 7199.61 7772.97 7643.62
		Water requirements
	Consumption (km³/year)	14107.86 14143.04 14143.04 14154.36 14154.36
	Withdrawal (km³/year)	61351.69 62191.48 62191.48 62228.12 62228.12
	CO 2 emissions (MtCO 2 /year) 38.20	34.43	34.38	28.18	27.46

Lowest flow on seven consecutive days for 10-year return period.

Flow with 95% of permanence over a period.

This is expressed in m 3 of fresh feed per m 3 of catalyst per hour. The inverse of LHSV is generally called residence time(STANISLAU et al, 2010).

Acknowledgements

We would like to acknowledge CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) and CNPq (National Council for Scientific and Technological Development) for financial support in the form of scholarships. We also thank IFP Energies nouvelles.


Luis Almeida
email: luis.almeida@sorbonne-universite.fr

Jorge Estrada

Nicolas Vauchelet
email: vauchelet@math.univ-paris13.fr

The sterile insect technique used as a barrier control against reinfestation

The sterile insect technique consists in massive release of sterilized males in the aim to reduce the size of mosquitoes population or even eradicate it. In this work, we investigate the feasability of using the sterile insect technique as a barrier against reinvasion. More precisely, we provide some numerical simulations and mathematical results showing that performing the sterile insect technique on a band large enough may stop reinvasion.

Introduction

Due to the number of diseases that they transmit, mosquitoes are considered as one of the most dangerous animal species for humans. According to the World Health Organization [1], vectorborne diseases account for more than 17% of all infectious diseases, causing more than 700 000 deaths annually. More than 3.9 billion people in over 128 countries are at risk of contracting dengue, with 96 millions cases estimated per year.

Several of the most prominent mosquito-transmitted diseases do not have a widely available vaccine; dengue, for instance, has only one licensed vaccine, Dengvaxia, with limited application [START_REF] Redoni | Dengue: Status of current and under-development vaccines[END_REF]. Zika and chikungunya do not have licensed vaccines, although several candidates are on diverse stages of development [23] [START_REF] Shan | Recent Progress in Vaccine Development Against Chikungunya Virus[END_REF]. Hence, at this stage, a natural strategy to control them is to act directly on the mosquito population. For this purpose, several strategies have been developed and experimented. Some techniques aim at replacing the existing population of mosquitoes by a population unable to propagate the pathogens, using the Wolbachia bacteria [START_REF] Hoffmann | Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission[END_REF]. Other techniques aim at reducing the size of the mosquito population : The Sterile Insect Technique (SIT) [START_REF] Dyck | Sterile Insect Technique Principles and Practice in Area-Wide Integrated Pest Management[END_REF][START_REF] Alphey | Sterileinsect methods for control of mosquito-borne diseases: an analysis[END_REF], the release of insects carrying a dominant lethal (RIDL) [START_REF] Thomas | Insect population control using adominant, repressible, lethal genetic system[END_REF][START_REF] Heinrich | A repressible female-specific lethal genetic system for making trans-genic insect strains suitable for a sterile-release program[END_REF][START_REF] Fu | Female-specific flightless phenotype for mosquito control[END_REF], the driving of anti-pathogen genes into natural populations [START_REF] Gould | A Killer Rescue system for self-limiting gene drive of anti-pathogen constructs[END_REF][START_REF] Marshall | Semele: a killer-male, rescuefemale system for suppression and replacement of insect disease vector populations[END_REF][START_REF] Ward | Medea selfish genetic elements as tools for altering traits of wild populations: a theoretical analysis[END_REF]. Finally, other strategies combine both reduce and replace strategies [START_REF] Robert | A reduce and replace strategy forsuppressing vector-borne diseases: insights from a deterministic model[END_REF].

In this article, we focus on the SIT. This strategy has been introduced in the 50's by Raymond C. Bushland and Edward F. Knipling. It consists in using area-wide releases of sterile insects to 1 reduce reproduction in a field population of the same species. Indeed wild female insects of the pest population do not reproduce when they are inseminated by released, radiation-sterilized males. For mosquito populations, this technique and the closely related incompatible insect technique (IIT) has been successfully used to drastically reduce their size in some isolated regions (see e.g. [START_REF] Stoll | Suppression of an isolated popu-lation of the mosquito vector Aedes polynesiensis on the atoll of Tetiaroa, French Polynesia, by sustained release of Wolbachia-incompatible male mosquitoes[END_REF][START_REF] Zheng | Incompatible and sterile insect techniques combined eliminate mosquitoes[END_REF]). In order to predict the dynamics of mosquito populations, mathematical modeling is an important tool. In particular, there is a growing interest in the study of control strategies (see e.g. [START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF][START_REF] Li | Modelling releases of sterile mosquitoes with different strategies[END_REF][START_REF] Bliman | Implementation of Control Strategies for Sterile Insect Techniques[END_REF][START_REF] Almeida | Control strategies on mosquitos population for the fight against arboviruses[END_REF] and references therein). in order to conduct rigorous studies, these works usually neglect spatial dependency.

In fact, only few articles propose to incorporate the spatial variable in their study of SIT. In [START_REF] Lewis | Wave-block in excitable media due to regions of depressed excitability[END_REF] the authors propose a simple scalar model to study the influence of the sterile insects density on the velocity of the spatial wave of spread of mosquitoes. The work [START_REF] Oléron Evans | A spatial model with pulsed releases to compare strategies for the sterile insect technique applied to the mosquito Aedes aegypti[END_REF] focuses on the influence of the release sites and the frequency of releases in the effectiveness of the sterile insect technique. In [START_REF] Seirin Lee | Optimal barrier zones for stopping the invasion of Aedes aegypti mosquitoes via transgenic or sterile insect techniques[END_REF][START_REF] Seirin Lee | Modelling Aedes aegypti mosquito control via transgenic and sterile insect techniques: Endemics and emerging outbreaks[END_REF], the authors conduct a numerical study on some mathematical models with spatial dependency to investigate the use of a barrier zone to prevent invasion of mosquitoes. However, up to our knowledge, there are no rigorous mathematical results on the existence of such barrier zone. In this paper, we conduct a similar study as in [START_REF] Seirin Lee | Optimal barrier zones for stopping the invasion of Aedes aegypti mosquitoes via transgenic or sterile insect techniques[END_REF] for another mathematical model which has been recently introduced in [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF]. Moreover, we propose a strategy to rigorously prove the existence of a barrier zone under some conditions on the parameters.

The outline of the paper is the following: In subsection 2.1 we introduce our dynamical system model and describe the variables and biological parameters, and also a simplified model that results from additional assumptions. We analyze the existence of positive equilibria and the stability of the mosquito-free equilibrium. In subsection 2.2 we introduce spatial models including diffusion. In section 3 we perform numerical simulations for the spatial models to observe the existence of wave-blocking for a large enough release of sterile males. In section 4 we prove the existence of a blocking phenomenon (some useful technical results are postponed to the appendix). We end this paper with a conclusion.

Mathematical model

Dynamical system

In the recent paper [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF], the following mathematical model governing the dynamics of mosquitoes has been proposed:

dE dt = b(1 - E K )F -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dF dt = rν E E(1 -e -β(M +γMs) ) M M + γM s -µ F F, dM s dt = u -µ s M s .
(2.1)

In this system, the mosquito population is divided into several compartments. The number of mosquitoes in the aquatic phase is denoted E; M and F denote respectively the number of adult wild males and adult females which have been fertilized; M s is the number of sterile mosquitoes which are released, the release function being denoted u. The fraction M M +γMs corresponds to the probability that a female mates with a wild male. Moreover, the term (1 -e -β(M +γMs) ) has been introduced to model a strong Allee effect, that is, the fact that when the number of males is close to zero, it will be very difficult for a female find a wild male to mate with.

Finally, we have the following parameters :

• b > 0 is the oviposition rate;

• µ E > 0, µ M > 0, µ F > 0, and µ s > 0 denote the death rates for the mosquitoes in the aquatic phase, for wild adult males, for adult females, and for sterile males, respectively;

• K is an environmental capacity for the aquatic phase, taking also into account the intraspecific competition;

• ν E > 0 is the rate of emergence;

• r ∈ (0, 1) is the probability that a female emerges and (1 -r) is the probability that a male emerges;

• u is a control function, corresponding to the number of sterile males which are released into the field.

Before introducing the spatial dependance in this system, we will first investigate the equilibria of (2.1) and their stability. For future use, we introduce the notation

g(F, M s ) = rν E bF b K F + µ E + ν E 1 -exp -β (1 -r)ν E bF µ M ( b K F + µ E + ν E ) + γM s (2.2) × (1 -r)ν E bF (1 -r)ν E bF + γµ M M s ( b K F + µ E + ν E ) -µ F F.
Then, we may use Lemma 3 in [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF] to prove the following result:

Proposition 2.1 Let N = brν E µ F (ν E +µ E ) and ψ = µ M
(1-r)ν E Kβ . Let us assume that N > max{1, 4ψ}. Let θ 0 ∈ (0, 1) be the unique solution of 1 -θ 0 = -4ψ N ln θ 0 and assume that max

θ∈[θ 0 ,1] -ln θ - 1 2ψ 1 -1 + 4ψ N ln(θ) 1 -θ > 0. (2.3)
Then, we have :

1. When u = 0. There exist two positive equilibria for system (2.1). They are given by

(E 1 , M 1 , F 1 , 0) and (E 2 , M 2 , F 2 , 0), with F 1 < F 2 and E i = bF i b K F i +ν E +µ E , M i = (1-r)ν E µ M E i , for i = 1, 2.
2. There exists a positive constant Ũ large enough, such that if u = U is a constant such that U > Ũ , then the unique equilibrium for system (2.1) is the mosquito-free equilibrium (0, 0, 0, U /µ s ), which is globally stable.

Proof. At equilibrium, we have

E = bF b K F + ν E + µ E , M = (1 -r)ν E µ M E.
Injecting into the equilibrium equation for F , we deduce that any equilibrium should satisfy g(F, M s ) = 0.

(2.4)

For u = 0, the only equilibrium of

dMs dt = u -µ s M s is M s = 0. Substituting M s = 0 in (2.
2), we obtain for 2.4

g(F, 0) = rν E bF b K F + µ E + ν E 1 -exp -β (1 -r)ν E bF µ M ( b K F + µ E + ν E ) -µ F F = 0.
F = 0 is always a solution. When F > 0 the equation g(F, 0) = 0 may be rewritten

g 1 (F ) := rν E b 1 -exp -β (1 -r)ν E bF µ M ( b K F + µ E + ν E ) = µ F ( b K F + µ E + ν E ).
The right hand side is an affine function. It is straightforward to verify that g 1 is an increasing concave function on (0, +∞) such that g 1 (0) = 0 and lim

F →+∞ g 1 (F ) = rν E b 1 -exp - β(1 -r)ν E K µ M .
Then, this concave function may intersect the affine function in 0,1 or 2 points. Notice that a necessary (but not sufficient) condition for intersection is that lim

F →+∞ g 1 (F ) > µ F (µ E + ν E ). The fact that condition (2.
3) guarantees existence of two positive roots is proved in [30, Lemma 3].

For the second part, for a given constant u = U , the only equilibrium of dMs dt = u -µ s M s is M s = U /µ s . Then, injecting this relation into the expression for g in (2.2), we deduce that for any F ≥ 0, we have lim U →+∞ g(F, U µs ) = -µ F F . For any α ∈ (0, µ F ), we may find U large enough such that for any F ≥ 0, g(F, U µs ) ≤ -αF . Hence, the only solution of the equation g(F, U µs ) = 0 is F = 0, which implies that the only equilibrium for system (2.1) is the mosquito-free equilibrium (0, 0, 0, U /µ s ).

Moreover, the Jacobian at the mosquito-free equilibrium is given by the following matrix

J(0, 0, 0, U /µ s ) =     -(ν E + µ E ) 0 b 0 (1 -r)ν E -µ M 0 0 0 0 -µ F 0 0 0 0 -µ s     .
We compute easily its eigenvalues which are given by {-ν E -µ E , -µ M , -µ F , -µ s }; they are all negative. Therefore the mosquito-free equilibrium is globally stable.

Remark 2.2 Although the proof does not give us an analytic formula for the equilibria, remark 3.2 in [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF] gives us a useful approximation of M 2 when β is not too small, as

M 2 ≈ 1 ψβ (1 -1 N ).

Spatial model

In order to model the spatial dynamics, we consider that adult mosquitoes diffuse according to a random walk. It is classical to model this active motion by adding a diffusion operator in the compartment of females and males. We denote by x the spatial variable. In order to simplify the approach, we only consider the one dimensional case, x ∈ R. Then, all unknown functions depend now on time t > 0 and position x ∈ R. In this setting, model (2.1) becomes

∂ t E = b(1 - E K )F -(ν E + µ E )E, (2.5a)
∂ t M -D u ∂ xx M = (1 -r)ν E E -µ M M, (2.5b)
∂ t F -D u ∂ xx F = rν E E(1 -e -β(M +γMs) ) M M + γM s -µ F F, (2.5c
)

∂ t M s -D u ∂ xx M s = u -µ s M s . (2.5d)
In this model D u is a given diffusion coefficient, which is assumed to be the same for both male and female mosquitoes. This system is complemented with some nonnegative initial data denoted E 0 , M 0 , F 0 , and M 0 s = 0. Then, there exists a nonnegative solution of (2.5). From now on, we will always assume that the parameters of the model are such that there are two positive equilibria for system (2.1) (see Proposition 2.1).

In the case when no control technique is implemented, i.e. M s = 0, we expect that this model to predict the invasion of the whole domain by the wild mosquitoes. In order to illustrate this phenomenon, we perform some numerical simulations of system (2.5). This system is discretized thanks to a semi-implicit finite difference scheme on an uniform mesh. We use the numerical values in Table (1), which are taken from [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF]. When M s = 0, the numerical results are shown in Figure 1. We observe that there is a wave of mosquitoes that invades the whole domain. Hence, in the absence of control and in a homogeneous environment, it is expected that the solutions to (2.5) invade freely the spatial domain.

If we fix a constant value of the number of sterile males in the whole domain, then the invasion is expected to slow down. This observation has already been studied in [START_REF] Lewis | Waves of extinction from sterile insect release[END_REF]. Numerical illustrations of this phenomenon are provided in Figure 2. In this simulation, the number of sterile males is fixed to be M s = 5000 on the whole domain. This figure should be compared to Figure 1 where M s = 0 : we observe that the invasion is slowed down.

In this paper, we want to investigate the possibility of blocking the propagation of the spreading of mosquitoes by releasing sterile mosquitoes on a band of width L. Can the sterile insect technique be used to act as a barrier to avoid re-invasion of wild mosquitoes in a mosquito-free region ? In order to answer this question, we first perform some numerical simulations in the next section.

Numerical simulations

We choose u(t, x) = U 1 [0,L] (x), where U is a given positive constant. We propose some numerical simulations. As above, we implement a finite difference scheme on an uniform mesh. The values of the numerical parameters are taken from [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF] and are given in Table 1.

Parameter β b r µ E ν E µ F µ M γ µ s D u Value
10 -2 10 0.49 0.03 0.05 0.04 0.1 1 0.12 0.0125

Table 1: Table of the numerical values used for the numerical simulations. These values are taken from [START_REF] Strugarek | On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations[END_REF] We show in Figures 3 the dynamics in time and space of the fertilized female population F for model (2.5). The domain where the release of sterile males is performed is of width L = 5 km. The release intensity is U = 10 000 km -1 day -1 (left), U = 20 000 km -1 day -1 (center), U = 30 000 km -1 day -1 (right). These results illustrate the influence of the intensity of the release U and of the width of the band on the efficiency of the blocking. We first notice that it seems that for U or L large enough, the mosquito wave is not able to pass through the release zone. On the contrary, if U and L are not large enough, the wave is only delayed by the release in the barrier zone.

It is also interesting to observe that when β → +∞, there is no blocking as it is illustrated in Figure 4 for model (2.5). In order to justify this observation, we observe that, when β → +∞ the mosquito steady state equilibrium is unstable for model (2.1). Indeed, when M s = 0 in system (2.1), after passing to the limit when β → +∞, the Jacobian matrix of the resulting system at the point Figure 3: Numerical simulations of system (2.5) with L = 5 km and U = 10 000 km -1 day -1 (left), U = 20 000 km -1 day -1 (center), and U = 30 000 km -1 day -1 (right). 1, while on the right we have the solution of the limit equation when β → +∞.

(0, 0, 0) is given by 

 -ν E -µ E 0 b (1 -r)ν E -µ M 0 rν E 0 -µ F  
Clearly -µ M is one of the eigenvalues. For brν E > µ F (ν E + µ E ), there is a positive real eigenvalue and therefore, the steady state would be unstable.

Mathematical approach

These numerical simulations seem to show that it is possible to block the spreading by releasing enough sterile males on a sufficiently wide domain. However, to be sure that it is not a numerical artifact and that the propagation is really blocked, we have to prove rigorous mathematical results.

The study of blocking by local action has been done by several authors with applications for instance in biology or criminal studies [START_REF] Pauwelussen | One way traffic of pulses in a neuron[END_REF][START_REF] Lewis | Wave-block in excitable media due to regions of depressed excitability[END_REF][START_REF] Chapuisat | Asymptotic profiles for a traveling front solution of a biological equation[END_REF][START_REF] Berestycki | Traveling wave solutions in a reaction-diffusion model for criminal activity[END_REF][START_REF] Nadin | Hindrances to bistable front propagation: application to Wolbachia invasion[END_REF][START_REF] Eberle | Front blocking versus propagation in the presence of a drift term in the direction of propagation[END_REF]. In this section, we will always consider that

u(t, x) = U 1 [0,L] (x),
and we assume that (2.3) is satisfied such that when u = 0 there are two stable steady states. We justify rigorously that for L > 0 and U large enough, we may block the propagation. We first provide some useful bounds and comparisons on the variables of our system. Then, we construct a stationary supersolution for (2.5) that vanishes at +∞. When it exists, we shall call this stationary supersolution a barrier for (2.5), since such a solution acts as a barrier to block the propagation.

Estimates

We first recall that, since initial data is nonnegative, the solution (E, F, M, M s ) of (2.5) is also nonnegative. Moreover, we have the following estimates : Lemma 4.1 Let us assume that E 0 ≤ K. Then, for any t ≥ 0, we have E ≤ K. Moreover, M and F are uniformly bounded.

Proof. The inequality is assumed to be true for t = 0. If there exists some t 0 > 0 and x ∈ R such that E(t 0 , x) ≥ K, then we deduce from the first equation in (2.5) that d dt E(t 0 , x) < 0, hence E(•, x) is decreasing in the vicinity of t 0 . Therefore, if this inequality is true initially, it should be true for any larger time. It follows then from a standard application of the maximum principle that M and F are also uniformly bounded.

Lemma 4.2 Let u(t, x) = U 1 [0,L] (x).
Then, the solution M s (t, x) to (2.5d) with initial data M s (t = 0) = 0 converges uniformly on x, when t → +∞, to the solution M s of the stationary equation

-D u ∂ xx M s = U 1 [0,L] -µ s M s , on R. (4.6)
Proof. We have

∂ t (M s -M s ) -∂ xx (M s -M s ) = -µ s (M s -M s ), M s (t = 0) = 0.
Then, using the heat kernel, we deduce

M s (t, x) = M s (x) - R 1 √ 4πt e -y 2 4t -µst M s (x -y) dy.
The stationary solution M s has the expression

M s (x) = U 2 √ D u µ s L 0 e -µs
Du |x-y| dy.

We can see from (4.1) that M s is symmetric with respect to x = L 2 and decreasing on ( L 2 , +∞). Therefore M s has a global maximum at x = L 2 , which we compute to be

M s ( L 2 ) = U µs 1 -e -µs Du L 2
, and we have

0 ≤ M s (x) -M s (t, x) = R 1 √ 4πt e -y 2 4t -µst M s (x -y) dy ≤ e -µst M s ( L 2 )
where we use the fact that the integral of the heat kernel equals 1. Therefore, M s (t, x) converges to M s exponentially in t and uniformly in x. Proof.

Using the symmetry of M s ,w e have that for any x ∈ [0,

L * ] M s (x) ≥ M s (0) = M s (L) if L > L * , and M s (x) ≥ M s (L * ) if L ≤ L * . If L > L * , we see that M s (0) = U 2µ s (1 -e -µs Du L ),
and, if L ≤ L * , we have that

M s (L * ) = U 2µ s e -µs Du L * (e µs Du L -1).
An immediate consequence of this result is the following :

Corollary 4.4 Let ε > 0 and L * > 0. Then, there exists U large enough such that, for any bounded nonnegative function M and for any x ∈ [0, L * ], we have

(1 -e -β(M (x)+γM s(x)) ) M (x) M (x) + γM s (x) ≤ ε,
where M s is the solution of (4.6).

Comparison principle

Let M s ≥ 0 be a given nonnegative function. Let us consider the stationary problem associated to (2.5) :

0 = b(1 - E K )F -(ν E + µ E )E (4.7a) -D u M = (1 -r)ν E E -µ M M (4.7b) -D u F = rν E E(1 -e -β(M +γM s) ) M M + γM s -µ F F . (4.7c)
We are particularly interested by stationary solutions which link the two stable equilibria (0, 0, 0) and (E 2 , M 2 , F 2 ) defined in Proposition 2.1. Therefore, we complement this system with conditions at infinity :

(E, M , F )(-∞) = (E 2 , M 2 , F 2 )
and (E, M , F )(+∞) = (0, 0, 0).

Definition 4.5

We call barrier any super-solution of (4.7) such that

(E, M , F )(-∞) = (E 2 , M 2 , F 2 )
and (E, M , F )(+∞) = (0, 0, 0).

We recall that super-solutions are obtained by replacing the equality signs in (4.7) by the inequalities ≥.

The following result explains why we call barrier such a solution : If the initial data is below this barrier, the solution stays below the barrier and hence the propagation (invasion of the whole domain by the wild population) will be blocked. Let us assume that there exists a barrier (E, M , F ) in the sense of definition 4.5.

If we assume that E 0 ≤ E, M 0 ≤ M , F 0 ≤ F . Then, the solution (E, M, F ) of (2.5a)-(2.5c) is such that for any t ≥ 0 and any x ∈ R, we have E(t, x) ≤ E, M (t, x) ≤ M and F (t, x) ≤ F .

Proof. From Lemma 4.1, we may consider the system (2.5a)-(2.5c) on the set [0

, K] × R + × R + . Let us denote g(M, M s ) := M M + γM s (1 -e -β(M +γMs) ).
We verify easily that, on this set, the map g is nonincreasing with respect to M s and nondecreasing with respect to M . Therefore, the system is monotonous and we deduce directly the result from standard arguments.

For the sake of completeness, we explain briefly these arguments. Let us denote δE = E -E, δM = M -M , and δF = F -F . Subtracting and (4.7) from (2.5), we deduce

∂ t δE ≤ b(1 - E K )δF -( bF K + ν E + µ E )δE ∂ t δM -D u ∂ xx δM ≤ (1 -r)ν E δE -µ M δM ∂ t δF -D u ∂ xx δF ≤ rν E δEg(M, M s ) + rν E E(g(M, M s ) -g(M , M s )) -µ F δF,
where we use the fact that g(M , M s ) ≤ g(M , M s ). Denoting by v + := max{v, 0} the positive part of any real v, we multiply the first equation by (δE) + , the second by (δM ) + , and the third by (δF ) + .

Integrating over R, we get 1 2

d dt R (δE) 2 + dx ≤ R b 1 - E K (δF ) + (δE) + dx - R bF K + ν E + µ E (δE) 2 + dx 1 2 d dt R (δM ) 2 + dx + D u R |∂ x (δM ) + | 2 dx ≤ R (1 -r)ν E (δE) + (δM ) + dx - R µ M (δM ) 2 + dx 1 2 d dt R (δF ) 2 + dx + D u R |∂ x (δF ) + | 2 dx ≤ R rν E (δE) + (δF ) + dx - R µ F (δF ) 2 + dx + R rν E E g(M, M s ) -g(M , M s ) (δF ) + dx,
where we use also the fact that g(M, M s ) ≤ 1. Since M → g(M, M s ) is increasing and uniformly Lipschitz-continuous, we have for some nonnegative constant

C g(M, M s ) -g(M , M s ) ≤ C(δM ) + .
Hence, adding these latter inequalities, we deduce that for some nonnegative constant C, we have

d dt R (δE) 2 + + (δF ) 2 + + (δM ) 2 + dx ≤ C R (δE) 2 + + (δF ) 2 + + (δM ) 2 + dx.
We conclude by using a Gronwall lemma, noticing that from our assumption on the initial data we have that at time t = 0, δE(t = 0) = 0, δM (t = 0) = 0, δF (t = 0) = 0.

Existence of a barrier

Let L * > 0 and α ∈ (0, µ F ), we deduce from Corollary 4.4 that for any L > 0, there exists U large enough such that for all x ∈ [0, L * ] the solution M s of (4.6) verifies, for any nonnegative and bounded functions M and F ,

rν E F (x) b K F (x) + µ E + ν E (1 -e -β(M (x)+γM s(x)) ) M (x) M (x) + γM s (x) ≤ (µ F -α)F (x).
It invites us to introduce the stationary problem

E = bF b K F + µ E + ν E (4.9a) -D u M = (1 -r)ν E E -µ M M (4.9b) -D u F = -αF 1 [0,L * ] + rν E E(1 -e -βM ) -µ F F 1 [0,L * ] c . (4.9c)
Any super-solution of (4.9) which satisfies the boundary conditions

(E, M , F )(-∞) = (E 2 , M 2 , F 2 )
and (E, M , F )(+∞) = (0, 0, 0) is a barrier, as defined in Definition 4.5. Then, in order to prove existence of a barrier, it suffices to prove that there exist super-solutions of (4.9), i.e. (M , F ) such that

-D u M ≥ (1 -r)ν E bF b K F + µ E + ν E -µ M M , on R \ {0, L}, (4.10a)
-D u F ≥ -αF 1 [0,L] + rν E bF b K F + µ E + ν E (1 -e -βM ) -µ F F 1 [0,L] c , on R \ {0, L}, (4.10b) M (ξ -) ≥ M (ξ + ), F (ξ -) ≥ F (ξ + ), for ξ ∈ {0, L}, (4.10c)
(M , F )(-∞) = (M 2 , F 2 ), (M , F )(+∞) = (0, 0). (4.10d)

We recall that (M 2 , F 2 ) is the larger mosquito equilibrium defined in (2.1). The main result of this work is the following :

Theorem 4.7 Under the same assumptions as in Proposition 2.1, for any L * > 0, there exists U large enough, such that there exists a solution (M , F ) of (4.10).

Hence, there exists a barrier for system (2.5).

Proof. Let α ∈ (0, µ F ) be fixed. We take U and L large enough such that Lemma 4.6 holds. The last point of the Theorem is a consequence of the comparison principle. Therefore, we just have to construct a super-solution of (4.9) which verifies the boundary conditions. We construct this solution piecewise.

First, on (-∞, 0), we notice that system (4.9) is the stationary system for (2.5) when M s = 0. From Proposition 2.1, we deduce that we may take M , F on (-∞, 0) as the larger values of the equilibria in 2.1, that is, M = M 2 and F = F 2 .

On (0, L), we introduce the notation τ := ν E b µ E +ν E and take (M (L), F (L)) such that

0 < M (L) < - 1 β ln 1 - µ F 2rτ , 0 < F (L) < µ M (1 -r)τ M (L). (4.11)
Theorem 4.7 establishes that it is possible to block a propagating front of invading mosquitoes by performing a sterile insect technique on a domain wide enough with releases of a sufficiently large amount of sterile males. However, the mathematical proof presented above provides only a sufficient condition to block the front and it is difficult to quantify the width of the domain and the number of sterile males. This is due to the fact that here we have to deal with a system of differential equations, whereas for a scalar equation it is possible to obtain a better characterization of the barrier.

Conclusion

In this work we have proved that mosquito invasion may be blocked by releasing a sufficient amount of sterile males over a sufficiently wide band. The main motivation of this result is to use the sterile insect technique to build a sanitary cordon to protect a certain region (e.g. an urban area) from wild mosquitoes living in an exterior region (a reservoir area where they are abundant like, for instance, a forest) and which are vectors of a disease (like dengue, zika or any other mosquito borne disease). We remark that the mathematical result is not exclusive to mosquitoes. This technique of building a barrier to avoid invasion may be applied to other invading insect species (for instance, other disease vectors or agricultural pests). It is however important that the population dynamics be characterized by an Allee effect such that the extinction steady state is stable.

Since we are looking for a solution v which is decreasing, we may invert it into a function x(v). Let us denote w(v) = -v (x(v)). Then, we have

w (v) = (µ -ψ(x(v)))v w(v) , w(v 0 ) = -v 0 .
The question is to know whether there exists v 0 < 0 such that the solution of this equation is defined on (0, v 0 ) and such that w(0) = 0. There are two possibilities: either w vanishes on v c ∈ (0, v 0 ) and then the solution is defined only on (v c , v 0 ) and we say that this solution is type I, or w does not vanish on (0, v 0 ) and then we set v c = 0 and we say that this solution is type II. We are looking for a solution which is between type I and type II.

Clearly, if v 0 < v 0 , then, from the Cauchy-Lipschitz theorem the corresponding solutions verify w < w. Therefore, the map v 0 → v c is nonincreasing and continuous. Moreover, we have seen at the beginning of this proof that if v 0 > -µ 2 v 0 the solution is type I and if v 0 < -√ µv 0 the solution is type II. Thus, by continuity there exists v 0 such that v c = 0 and w(0) = 0. Moreover, we have also v 0 ∈ [-√ µv 0 , -µ 2 v 0 ] and with (A.5) we deduce the estimate v(x) ≤ 1 2 (v 0 + v 0 √ µ )e -√ µx .

In the same spirit we have the following result:

Lemma A.3 Let v be a nonincreasing and nonnegative function on [0, ∞). Let u 0 > 0. We assume that αv(0) < βu 0 . We consider the Cauchy problem on (0, +∞) where we use the fact that v is nonnegative. Hence the solution is type I if u 0 < -√ βu 0 . Moreover, using the fact that v(x) ≤ v(0), we have u(x) ≥ u 0 cosh( βx) + Therefore, the solution is type II if u 0 > 1 √ β (αv(0) -βu 0 ). By a continuity and monotonicity argument, as in the proof of Lemma A.2, there exists u 0 ∈ [-√ βu 0 , 1 √ β (αv(0) -βu 0 )] such that the solution of (A.6) is such that lim x→+∞ u(x) = 0. Moreover, we have from the first estimate in this proof that u(x) ≤ 1 2 (u 0 +

-u = αv -βu, u (0)
u 0 √ β )e √ βx .

Figure 1 :

 1 Figure 1: Numerical simulation of the time and space dynamics of (F, M ) solving (2.5) when M s = 0. Left : Profiles of solutions at three different times with same initial data; solutions are plotted at time t = 100, t = 200, t = 300. Right : Dynamics in time and space of the fertilized female density for system (2.5).

Figure 2 :

 2 Figure 2: Numerical simulation of the time and space dynamics of (F, M ) solving (2.5) when M s = 5000. Left : Profiles of solutions at three different times with same initial data; solutions are plotted at time t = 100, t = 200, t = 300. Right : Dynamics in time and space of the females density for system (2.5). Comparing this figure with Fig. 1, we observe that the invasion is slowing down.

Figure 4 :

 4 Figure 4: Numerical simulations for the model (2.5) with β → +∞, L = 10 km and U = 40 000 km -1 day -1 . In the figure on the left β = 10 -2 as in Table1, while on the right we have the solution of the limit equation when β → +∞.

Lemma 4 . 3

 43 Let L > 0. Then, for any L * > 0, the solution M s of (4.6) verifies, for any x ∈ [0, L * ]

Lemma 4 . 6

 46 Let M s ≥ M s . Let us assume that the initial data are nonnegative and such that E 0 ≤ K.(4.8)

  = u 0 , u (0) = u 0 . (A.6)Then, there exists u 0 ∈ [-√ βu 0 , 1 √ β (αv(0) -βu 0 )] such that the solution of (A.6) is decreasing and tends to 0 at +∞. Moreover, u(x) ≤ Ce - √ βx for some nonnegative constant C.Proof. Following the idea of the proof of Lemma A.2, we are left to prove that there exist type I and type II solutions. From the Duhamel formula, we haveu(x) = u 0 cosh( βx) + u 0 √ β sinh( βx) -α √ β x 0 v(z) sinh( β(x -z)) dz ≤ u 0 cosh( βx) + u 0 √ β sinh( βx),

√ µ M Du ((1 -r)τ F (L) -µ M M (L)).It concludes the construction of (M , F ).

Acknowledgment

The authors acknowledge the unknown referee for his/her comments and suggestions which allowed to improve this paper.

Then, we solve (4.9) on (0, L) with boundary conditions (M (0), F (0)) = (M 2 , F 2 ) and (M (L), F (L)). From (4.9c), we deduce that

which is positive on (0, L). From (4.9b), we deduce using the comparison principle that M ≥ M > 0, where M is the solution of -D u M + µ M M = 0 with M (0) = M 2 and M (L) = M (L).

On (L, +∞), we first observe that for F ≥ 0, we have

prove that there exists a nonnegative solution on (L, +∞), that tends to 0 at +∞, of the system

with (M (L), F (L)) as in (4.11). Notice that from (4.11), we have that (A.3) holds with φ(u) = rτ (1 -e -βu ) (which is clearly increasing and Lipschitz continuous). Then, we may apply the result in Proposition A.1 in the Appendix : For F (L) small enough, there exists (M (L + ), F (L + )) such that the solution of the above Cauchy problem is nonincreasing, nonnegative and goes exponentially fast to 0 as x tends to +∞. Thus, our solution so far satisfies (4.10a) and (4.10b), and it only remains to prove that (4.10c) holds in order to conclude our proof.

From the expression of F on (0, L), we compute

We already know that F (0 -) = 0, and from Lemma A.2, we have

In the same manner M (0 -) = 0, and we have from Lemma A.

. By the comparison principle, M is bounded from above by the solution on (0, L) of

Then, we have M (0 + ) ≤ y (0) and M (L -) ≥ y (L) where, after some straightforward computations,

We can always find L large enough and (M (L), F (L)), with F (L) small, such that (4.11) holds true,

0) ≤ 0, and y (L) ≥

Appendix Existence of nonincreasing solution vanishing at infinity

In this appendix, we consider the second order differential system on (0, +∞)

We assume φ nonnegative, nondecreasing, and Lipschitz continuous , (A.2)

Our main result is the following Proposition A.1 Under assumptions (A.2) and (A.3), for v 0 small enough, there exist (u 0 , v 0 ) such that the solution (u, v) to (A.1) is such that u and v are nonincreasing and go exponentially fast to 0 as x → +∞.

Before proving this result, we will prove two technical lemmas which concern (A.1a) and (A.1b) independently.

Lemma A.2 Let µ > 0 and ψ be a nonincreasing and nonnegative function such that ψ(0) ≤ µ 2 , and let v 0 > 0. Let us consider the Cauchy problem on (0, +∞)

Then, there exists v 0 ∈ [-√ µv 0 , -µ 2 v 0 ] such that the solution of (A.4) is decreasing and tends to 0 at +∞. Moreover, v(x) ≤ Ce -√ µx for some nonnegative constant C.

Proof. From the Duhamel formula, we have that

as long as v ≥ 0. Hence, we deduce that if v 0 < -√ µv 0 , the right hand side vanishes which implies that v vanishes. From the assumptions on ψ, we also have that

Then, if v 0 > -µ 2 v 0 , we have v > 0 and v(x) goes to +∞ as x → +∞.

Proof of Proposition A.1 Let us consider a nonincreasing and nonnegative function u on [0, +∞) such that u(0) = u 0 . Then, from Lemma A.2, there exists v 0 such that the solution of

is nonincreasing, nonnegative, and decays exponentially fast to 0 at +∞. With such a function v, from Lemma A.3, there exists u 0 such that the solution of (A.6) is nonincreasing, nonnegative, and decays exponentially fast to 0 at +∞. We denote by F(u) this solution, which allows to define the map F. Let us denote A the closed subset of functions in H 1 (0, +∞) which are nonnegative and nonincreasing and take the value u 0 at 0. Clearly F maps A into itself. Let us prove that F is a contraction on A for v 0 small enough.

Let u 1 and u 2 be in A. Then, by definition of F, we have

On the one hand, multiplying (A.7) by (F(u 1 ) -F(u 2 )) and integrating, we obtain, after an integration by parts

We deduce from a Cauchy-Schwarz inequality

On the other hand, multiplying (A.8) by (v 1 -v 2 ) and integrating we get, in a similar way,

where we use (A.3), the Lipschitz continuity of φ (see (A.2)), and the fact that v 2 is decreasing. Applying again a Cauchy-Schwarz inequality, we obtain

With (A.9), we conclude that

Hence, F is a contraction for v 0 small enough. Thus, there exists a unique fixed point, which is a solution of our problem.