Nicolas Lamaison

Simon Colle

Mathieu Vallée

Roland Bavière

Storage Influence in a Combined Biomass and Power-to-Heat District Heating Production Plant

Keywords: District Hea ng, Power-to-Heat, Biomass, Storage, MILP, Model Predic ve Control, Numerical Simulator

The increasing penetra on of wind and solar electricity becomes challenging for grid operators. Interconnec ng electrical and thermal networks through the Power-to-Heat concept brings flexibility to the electrical grid while supplying a significant renewable source to District Hea ng (DH). In the present paper, we study a DH produc on plant composed of a biomass generator, a heat-pump and a heat storage in the French energe c context. We assess the techno-economic performances of this system using Mixed Integer Linear Programming (MILP). A mul -objec ve parametric op miza on method is applied to size the system using the available quan ty of biomass, the maximum CO2 content and minimum renewable energy ra o (REnR) of the heat produc on as ε-constraints. Our analysis shows that without strong constraints, heat pump and daily storage are used. For a limited amount of biomass available, we also show that inves ng in an inter-seasonal storage is necessary to reach high REnR. For comparisons, this energy system is also assessed with the Danish and German electric mix. We then verify the sizing stage results on the opera onal performances of a nonlinear numerical simulator. With that methodology, it is possible to evaluate the impact of the MILP modelling level of detail on the obtained results. An error of 5.1% on the produc on trajectories is here obtained for a given system.

Nomenclature

Greek Letters

Introduction

In France, similarly to other European countries, 35% of the final energy consump on is devoted to space hea ng and domes c hot water prepara on (DHW), which amounts to 665 TWh per year [START_REF]An EU Strategy on Hea ng and Cooling[END_REF]. The French energy planning of 2016 [START_REF] Francaise | Programma on Pluriannuelle de l'Energie[END_REF] has iden fied District Hea ng Networks (DHN) as a solu on to reduce the use of fossil energy to supply heat demands. Indeed, due to their ability to massively distribute renewable and recovery energies (R&R), DHN are expected to deliver 5 mes more R&R in 2030, to reach about 40TWh.

Similarly, the Heat Roadmap France [START_REF] Paardekooper | Heat Roadmap France: Quan fying the Impact of Low-Carbon Hea ng and Cooling Roadmaps[END_REF] recommends that DHN should cover up to 25% of the heat demand by 2050 while its current share is only 6%.

In parallel, the concept of 4 th Genera on District Hea ng [START_REF] Lund | 4th Genera on District Hea ng (4GDH)[END_REF][START_REF] Lund | The status of 4th genera on district hea ng: Research and results[END_REF] emphasizes the need to interconnect DHNs to the power grid, as a way to provide the flexibility required for a broader integra on of intermiGent renewable energies [START_REF] Vandermeulen | Controlling district hea ng and cooling networks to unlock flexibility: A review[END_REF]. In France, while combined Heat and Power (CHP) is limited by the rather low price of nuclearbased electricity, this integra on will then be first driven by Power-to-Heat (PtoH) [START_REF]Étude de valorisa on du stockage thermique et du power-to-heat[END_REF]. Indeed, combining PtoH with thermal storage provides a cheap and efficient form of storing excess renewable power [START_REF] Lund | Energy storage and smart energy systems[END_REF].

At the same me, and contrarily to other EU countries, biomass is expected to play a major role and reach 50% in the French DHN mix by 2030 [START_REF] Sncu | Enquête annuelle sur les réseaux de chaleur-res tu on des sta s ques 2014[END_REF]. It is worth men oning that the number of biomass-based DHN in France has significantly grown over the last 10 years, mainly because of financial incen ves. It is now es mated that about 500 DHN [START_REF] Amorce | Les réseaux de Chaleur au bois (rapport 2016 -données 2014)[END_REF] are using biomass in their energy mix. However, biomass must always be considered as a limited resource, unevenly distributed and affected by transporta on constraints [START_REF] Ericsson | Assessment of the poten al biomass supply in Europe using a resourcefocused approach[END_REF]. Moreover, its renewable nature depends on its usage rate.

In this paper, we study the op mal sizing of a DHN produc on plant combining PtoH, biomass and storage in the French context. In par cular, we consider the influence of several parameters on the need for different thermal energy storage sizes, from small water tanks to large inter-seasonal storages. We consider only waterbased storage at the produc on side, although other technologies as well as storing heat at building level or in the network itself could be considered at later stages [START_REF] Leśko | Opera onal op miza on in district hea ng systems with the use of thermal energy storage[END_REF].

According to the background and literature review, the paper is then structured as follows. In Sec on 2, we provide an overview of the state-of-the-art on similar systems, as well as on the typical sizing methodologies used for such studies. Sec on 3 presents the studied produc on plant, the methodology adopted in this study, which combines an op mal sizing stage with a sizing valida on stage. In Sec on 4, we detail the Mixed-Integer Linear Programming (MILP) model used for op mal sizing defini on, its implementa on and the results obtained under various annual constraints. Sec on 5 aims at valida ng the designed system by assessing its opera onal performance through simula on with a numerical simulator of the system in order to calculate the possible devia ons from the op mal sizing calcula on. In Sec on 6, we conclude and propose some perspec ves.

State-of-the-art

2.1.Smart energy systems combining biomass, power-to-heat and thermal storage capacities

As pointed out by Kwon and Ostergaard [START_REF] Koch | Intelligent Controlling of Power Driven Solid Biomass CHP Plants in Flexible District Hea ng With a Seasonal Heat Storage And a Power-To-Heat Component[END_REF], it seems rather evident that the biomass resource will be subjected to a sever supply stress in the coming years with for example beGer purposes elsewhere in the energy systems. Thus, methodologies that cannot account for the limited availability of this resource are inappropriate in prac ce.

The on-going work of Koch and Gaderer [START_REF] Koch | Intelligent Controlling of Power Driven Solid Biomass CHP Plants in Flexible District Hea ng With a Seasonal Heat Storage And a Power-To-Heat Component[END_REF] addresses the intelligent controlling of a system combining a soliddriven biomass cogenera on heat plant (CHP) and PtoH with heat storage capaci es enables a DHN to contribute to the electricity and balancing market, thereby increasing the overall system efficiency. In their recent work, a conceptual approach for studying the flexibility-oriented sizing of such system is described. Their methodology relies on a set of simulated scenarios ran on a numerical simulator but does not account for annual constraints limi ng.

The work of Ostergaard et al. [START_REF] Østergaard | Business and socioeconomic assessment of introducing heat pumps with heat storage in small-scale district hea ng systems[END_REF] inves gates the transi on from biomass-based DHN to heat pump and storage based DHN from a business and socio-economic point of view. It is highlighted that even though the heat pump and storage system is interes ng especially thanks to the flexibility it offers, it seems rather limited from an economical point of view. As pointed out in next sec on, it seems that the rule-based opera on they used to perform their analysis is not harnessing the full poten al of such energy system, thus reducing its economic benefit.

In the holis c approach of DHN produc on plant of Dahl et al. [START_REF] Dahl | Cost sensi vity of op mal sector-coupled district hea ng produc on systems[END_REF], applied in the Danish energe c context, it is shown that when the fossil-based resources are limited, inves ng in large-scale heat pumps and heat storages is desirable, phasing-out the use of CHPs, with only a slight increase of system cost.

As no ced previously, while biomass boiler is generally used as CHP throughout European DHNs, biomass boilers are seldom used as CHP in the French energe c context because of the largely nuclear-based electrical produc on. Even though there is a high poten al to combine the numerous available biomass boilers with PtoH and storage capaci es, no studies dealing with such a system in the French context could be found in the open literature.

2.2.Optimal sizing methodologies

Concerning the sizing approach, tradi onal methods use either dura on curves and linear cost characteris cs as explained in Frederiksen and Werner [START_REF] Frederiksen | District Hea6ng and Cooling[END_REF] or dynamic simula on tools with rule-based control strategies as used by Ostergaard et al. [START_REF] Østergaard | Business and socioeconomic assessment of introducing heat pumps with heat storage in small-scale district hea ng systems[END_REF] for combined heat pump and storage opera on and Lund et al. [START_REF] Frederiksen | District Hea6ng and Cooling[END_REF] for distributed small-scale CHP plants. Both methods are not capable of op mally sizing produc on plants when resources exhibit highly variable costs, e.g. electricity, or when the plant has storage capaci es. Moreover, adding annual constraints, such as CO2 or renewable contents of the produc on or limited annual quan ty of a given resource, e.g. biomass, to these methods is complex and me consuming since itera ve methods have to be used. On the other hand, opera onal performances of complex energy systems can be significantly underes mated when rule-based control is used. This has been shown by the work of Giraud et al. [START_REF] Giraud | Op mal control of district hea ng systems using dynamic simula on and mixed integer linear programming[END_REF] in the field of DHN opera on and Dufo-Lopez et al. [START_REF] Dufo-López | Daily opera on op misa on of hybrid stand-alone system by model predic ve control considering ageing model[END_REF] in the field of hybrid photovoltaics-based system opera on. These elements explain why the results obtained using both aforemen oned methods can be said sub-op mal.

In comparison, since their introduc on by Grossmann [START_REF] Grossmann | Mixed-integer programming approach for the synthesis of integrated process flowsheets[END_REF] for the study of energy systems, Mixed Integer Linear Programming (MILP) methods have been used successfully in numerous energy studies, e.g. for CO2 network and power to gas [START_REF] Suciu | Energy integra on of CO2 networks and power to gas for emerging energy autonomous ci es in Europe[END_REF], process integra on [START_REF] Marechal | Targe ng the integra on of mul -period u lity systems for site scale process integra on[END_REF], refrigera on systems [START_REF] Shelton | Op mal synthesis of integrated refrigera on systems-I: Mixedinteger programming model[END_REF] and solar assisted heat pumping [START_REF] Wallerand | Op mal design of solar-assisted industrial processes considering heat pumping: Case study of a dairy[END_REF]. A detailed review of op miza on studies in energy systems dealing with electricity, hea ng and cooling by Ünal et al. [START_REF] Ünal | Op misa on studies on tri-genera on: a review: Op misa on studies on tri-genera on: a review[END_REF], highlights the extensive growing usage of MILP approach in the community. The advantages of such a method are the following: i) the close rela on of a MILP formula on to the physical equa ons, which facilitates the expression of realis c models, ii) the reasonable computa onal me, and iii) the guarantee of finding global op mality, which especially ensures the perfect reproducibility of results when performing comparisons.

The main drawback is the requirement of using only linear or quasi-linear models. When dealing with energy systems, nonlinear behaviors are plen ful, e.g. fluid transporta on and temperature dependent efficiencies, leading to strong approxima on of MILP approaches in some cases. The laGer is the reason why some authors have inves gated global op mizers such as gene c algorithms, e.g. for trigenera on systems [START_REF] Giraud | Op mal control of district hea ng systems using dynamic simula on and mixed integer linear programming[END_REF][START_REF] Dufo-López | Daily opera on op misa on of hybrid stand-alone system by model predic ve control considering ageing model[END_REF] or Mixed Integer Non Linear Programming (MINLP), e.g. for mul -period CHP plants [START_REF] Tveit | Mul -period MINLP model for op mising opera on and structural changes to CHP plants in district hea ng networks with long-term thermal storage[END_REF]. Both strategies tend to be computa onally costly and thus not appropriate in sizing stage. Moreover, in both strategies, the available algorithms are not robust enough and may converge to local, instead of global, op mum.

Instead of trying to overcome the main drawback of MILP programming, i.e. the linear approxima on of all the physical phenomena, the present study uses a detailed non-linear numerical simulator of the system to validate and/or update the results from the MILP-based sizing. That innova ve methodology allows to incrementally add complexity to the sizing MILP un l sa sfactory opera on of the detailed numerical simulator is achieved. We present the adopted methodology in the next sec on of the present paper.

Methodology

The present sec on starts with the presenta on of the studied heat produc on plant. An overview of the methodology used throughout this study is then addressed. The boundary condi ons and input data are presented in sequence. Finally, the performance indicators we used are listed.

3.1.System description

Following the assessment presented in sec on 1, the present paper studies the produc on plant presented in Figure 1. It comprises a biomass boiler, a heat pump as PtoH element and a thermal storage in parallel arrangement. It also comprises a back-up gas boiler in serial arrangement just before the exit of the plant so that the demand is always sa sfied even during peak or maintenance periods.

3.2.Overview

Figure 2 gives an overview of the methodology used in this work, which combines an op mal sizing stage with a sizing verifica on stage.

Star ng from given input data (investment costs, efficiencies, technical limita ons) and boundary condi ons (DHN load, variable costs), op mal sizing of the system is firstly performed with quasi-linear models (Stage 1 in Figure 2). For the laGer, the decision variables are gathered in two groups, i.e. the opera onal ones that require a value at each me step (e.g. Power level of each equipment) and the sizing ones that have a unique value (e.g. size of the storage). For this sizing stage, a single op misa on is performed with a fixed horizon of 1 year. The results are then extrapolated over the life me of the produc on plant " # , set in the present study at 20 years. The results of this op mal sizing stage are i) the sizing of the system and ii) the opera onal variables trajectories and indicators (presented in sec on 3.4).

Figure 2: Overview of the methodology used in the present study

Secondly, the sizing is verified in opera on (Stage 2 in Figure 2) with the same given input data and boundary condi ons. For that stage, the sizing parameters are fixed and used by the opera onal control module. The laGer uses a Model Predic ve Control (MPC) algorithm to calculate the opera onal variables trajectories on a receding horizon. Results of this calcula on are used as set points for a numerical simulator of the system, i.e. a detailed thermal-hydraulic and non-linear model of the produc on plant. The outputs of this numerical

Biomass Boiler

Water Storage Heat Pump

Back-Up

Biomass

Electricity

Gas

Heat

To DHN simulator, i.e. the opera onal variables trajectories and indicators, are compared to the ones obtained during the sizing stage to validate the sizing.

The same boundary condi ons are considered for the op mal sizing, the opera onal control and the numerical simulator. The laGer means for example that the load curve considered during the sizing stage is the same as the one considered for the numerical simulator, which is represen ng real life opera on. The laGer thus represents idealis c condi ons with no uncertain es in the predic on models. It is worth men oning that the present methodology is able to deal with uncertain es on the boundary condi ons.

While studies dealing with the sizing of energy systems generally consider only the sizing based on a single MILP op miza on using an horizon of 1 year (Stage 1 in Figure 2), the present work uses a non-linear numerical simulator to both verify the sizing and evaluate more realis cally the opera onal performances (Stage 2 in Figure 2).

3.3.Boundary conditions and input data

This sec on presents the main boundary condi ons and input data used for both the sizing and the valida on stages in sec ons 4 and 5. Part of this data is highly dependent on the country, so specific cases for France, Denmark and Germany have been considered (see sec on 3.3.2). Denmark and Germany energe c contexts are here used as comparisons with the French case.

Country-independent data

Electricity Price.

The electricity price considered is composed of a constant part, represen ng taxes and transporta on, and a variable part, i.e. the electricity exchange price. In France, the constant part is approximately equal to 45€/MWh. As for the variable price, it was extracted from the 2017 EPEX-Spot database [START_REF]EPEX-spot[END_REF], the most used electricity exchange market in Europe. It is worth men oning here that the same electricity costs were considered in all the energe c contexts. In reality, even though the day-ahead market prices are similar [START_REF]EPEX-spot[END_REF], that assump on is incorrect since transporta on and taxes are different. However, this assump on is acceptable here since i) the error introduced in the electricity cost is in the order of 15% and ii) the objec ve here is to compare the electric mix influence.

Network Load

We generated a representa ve load using a tool [START_REF] Provent | Smart Grid Solaire Thermique: Rapport d'étude sur les réseaux de chaleur existants et les réseaux adaptés aux Eco-quar ers[END_REF] which provide an hourly profile demand of a DHN accoun ng for space hea ng, DHW prepara on and network losses. We consider 3000 equivalent dwellings, for which the distributed monthly consump on is presented in Table 1. We used weather data for the year 2017 in Chambéry in France, and we reconstructed cold water temperature (for DHW) and ground temperature (for DHN losses) using the models of Burch and Christensen [START_REF] Burch | Towards development of an algorithm for mains water temperature[END_REF] and Kasuda and Archenbach [START_REF] Kasuda | Earth Temperature and Thermal Diffusivity at Selected Sta ons in the United States[END_REF] respec vely. A hea ng law was considered for the supply network temperature (65 to 95°C) and a return temperature of 45°C and 55°C was set respec vely during hea ng (01/10 to 15/05) and non-hea ng seasons. DHN overall losses were assumed to represent 10% of the DHN supplied energy and were recalculated on an hourly profile using the ground temperature. The peak network load >?%, (@ and total demand >?%,A (; A obtained are respec vely of 18 MW and 40.10 3 MWh. This leads to a (-) ra o, as defined by Equa on (1), of 0.25, where 8760 is the number of hours composing a year and which is a typical value for DHN. The network load was considered the same in all the energe c-context so that the total heat demand was preserved. Costs.

For the produc on plant hea ng equipment, Table 2 lists the values of the investment and produc on costs ( 7 and ; 5 ). For the heat pump, the produc on costs are obtained by dividing the sum of the variable and constant electricity prices by the heat pump coefficient of performance (

), considered to be 3 in the present work.

= 3 3 3 (2)
where the heat pump is defined as the ra o of thermal power provided to the DHN 3 3 over the electrical consump on of the heat pump 3 .

Concerning the storage,

7

L and 7 M respec vely the investment related to the energy and power sizing, are obtained from a study dedicated to PtoH and storage in France [START_REF]Étude de valorisa on du stockage thermique et du power-to-heat[END_REF] and from the work of Eames et al. [START_REF] Eames | The Future Role of Thermal Energy Storage in the UK Energy System: An Assessment of the Technical Feasibility and Factors Influencing Adop on[END_REF]. 7 L is calculated using a two-parts piecewise linear formula on, in which the boGom part relates to low capacity tank (3 €/kWhth) while the upper part relates to high capacity pit storages (0.8 €/kWhth). For Table 2 contains other data such as the ra o of minimum power to maximum power ! of all the hea ng equipment. The laGer aims at represen ng the technical limits of the different components which play a decisive role in the way the produc on plant is conducted (e.g. switch off of the biomass boiler when the load is below the minimum technical limit). The values were averaged from various manufacturers catalogues. The last two rows of Table 2 are related to the CO2 content and Renewable Energy Ra o (REnR) of each equipment. Although the study focuses on the French context, we also performed a comparison with the Danish and German contexts, especially regarding the CO2 content and REnR of electricity. Figure 3 presents the CO2 emissions in g/kWhel of the electricity produc on for the three countries. It was obtained using the annual electricity produc on planning from the European Network of Transmission System Operators for Electricity database of year 2017 [START_REF] Entso-E | [END_REF]. In order to convert these profiles to CO2 emissions in g/kWhth of the heat pump produc on, the values are divided by the of the heat pump.

Similarly, Figure 4 presents the REnR of the electricity produc on for the year 2017 in France, Germany and Denmark. In order to convert these profiles to the REnR of the heat pump produc on, Equa on (3) must be used.

& '& 3 3 ( ) = & '& 8 9, ( ) + -1 (3)
where & '& 3 3 ( ) is the thermal renewable energy ra o of the heat pump, & '& 8 9, ( ) is given by Figure 4.

The main conclusions to draw from the profiles of Figure 3 and Figure 4 are the following: i) France exhibits both a low CO2 content but also a low REnR because of its nuclear-based electrical produc on, ii) in Denmark, the very-intermiGent wind power backed up by coal-fired power plant leads to highly variable CO2 content and REnR both at much larger values than in France, and iii) Germany lies in between the French and Danish contexts, since its electricity produc on exhibits a significant share of photovoltaics also backed-up by coalfired power plants.

3.4.Indicators

The main indicator used here is the Levelized Cost of Energy (LCOE), calculated with Equa on (4) in €/kWhth. It represents the average minimum price at which the produced heat must be sold in order to break-even over the life me of the produc on plant " # .

= ∑ ∑ (5 ( ). ∆ % bcdeb fg (1 + ( ) Ahg = ijkd Afg (4)
Where (5 ( ) [kW] represents the DHN heat load at me step t (see sec on 3.3.1), l [h] and represents respec vely the me step and the number of me steps, ( [%] is the annualiza on factor (set to 5%), is the sum of the capital expenditure (CAPEX) and the annualized opera ng expenditure (OPEX) over the life me of the equipment " # .

is calculated using Equa on (5).

= 7 + m 1 (1 + ( ) @hg ( ; 5 + (; + 8( ) = ijkd @fg [START_REF] Lund | The status of 4th genera on district hea ng: Research and results[END_REF] Where 7 , 8( , ; 5 , and (; are respec vely the investment, maintenance, produc on and star ng costs. These terms are further detailed in the next sec on.

Two other indicators are used, i.e. the yearly thermal REnR and CO2 content of the heat produc on, as defined respec vely in Equa ons ( 6) and [START_REF]Étude de valorisa on du stockage thermique et du power-to-heat[END_REF]. It is worth men oning that from now on, all the decision variables for the op mal sizing or the opera onal control will be highlighted in bold font to ease the reading of the equa ons.

= ∑ ∑ . n S (X). l % dop fg % bcdeb fg ∑ (5 ( ). l % bcdeb fg (6)
& '& 3 = ∑ ∑ & '& 3 ( ). n S (X) % dop fg % bcdeb fg ∑ (5 ( ) % bcdeb fg (7)
Where the and & '& 3 represent the CO2 content and the thermal REnR of the produc on by the different equipment. Their values were listed in Table 2. Addi onally, n S ( ) [kW], represents the power value at me step t of equipment i. From now on, the REnR will always be the thermal one.

Optimal System Sizing

The present sec on deals with the sizing problem formula on and implementa on of the op mal sizing presented in Figure 2 of the previous sec on. The obtained results are then discussed.

Base MILP problem formulation

At the op mal sizing stage, we formulate a MILP problem, with an objec ve to find the vector of decision variables q = = rq s , … q u , q uvs , … , q T w solu on of the problem of Equa on (8), q being composed of con nuous (1, … , x) and integer (x + 1, … , ') variables.

min 9 = = . q |} ℎ • € ≤ . q ≤ & € ) < ≤ q ≤ ' < (8)
where c [n] is a vector of cost, [m x n], €[m] and & €[m] are respec vely the matrix and vectors of linear constraints, and ) < [n] and ' < [n] are respec vely the lower and upper bounds vector for the decision variables.

Objec ve Func on

The problem objec ve in Equa on ( 8) is the minimiza on of the cost func on which is calculated with Equa on (9) (as in Equa on ( 5)). Where 7 , 8( , ; 5 , and (; are respec vely the investment, maintenance, produc on and star ng costs. For confiden ality reasons, it is generally rather complicate to get correct values for the star ng costs (; . In order to prevent incorrect op miza on because of too high star ng costs, the laGer have been set to a negligible value and are only used to write the constraint on the minimum ON me of the component, as explained in sec on 4. [kW] represents the sizing values of respec vely the maximum energy content of the storage, the maximum power of the storage, the maximum power of equipment i, and ‰ S ( ) is an integer variable equal to 1 when the equipment just started and 0 otherwise.

Modelling assump ons

In the sizing MILP model, the following assump ons were made in order to have a tractable problem:

-The of the heat pump does not vary with the temperature or power levels of the heat pump; -Thermal stra fica on is not accounted for in the storage; -Resource costs are not affected by the opera on of the DHN. In reality, the electricity or biomass costs may be influenced by the large consump ons of many such produc on plants; -In the sizing stage, the heat loss of the storage are considered null because a constant value (required for the linear models) would give misleading results since the losses should be dependent on the size of the storage which is a decision variable. Heat loss are however accounted for in opera on in sec on 5.

Equality Constraints

Plant Energy Balance. As wriGen in Equa on [START_REF] Amorce | Les réseaux de Chaleur au bois (rapport 2016 -données 2014)[END_REF], the total DHN heat load must be met at each me step t.

m n S Š V‹OE Sfs (X) + n RSWN_ WX (X) = n N_ WX (X) + (5 ( ) (10)
where n RSWN_ WX (X) [kW] and n N_ WX (X) [kW] represent the discharging and charging power at me step t.

Storage Energy Balance. As wriGen in Equa on [START_REF] Ericsson | Assessment of the poten al biomass supply in Europe using a resourcefocused approach[END_REF], a second energy balance can be wriGen for the storage.

^WX (X) -^WX (X -s) l = n N_ WX (X) -n RSWN_ WX (X) - . ^WX (X) (11)
where ^WX (X) [kWhth] is the energy content of the storage at me t, and [h -1 ] is the overall heat loss coefficient. As explained beforehand, the laGer is set to null in the sizing phase but is used during the opera onal valida on in sec on 5.

Inequality Constraints

Equipment Power Capacity Bounds. The power supplied by each equipment is bounded by minimum and maximum values which are also decision variables, as shown in Equa on [START_REF] Leśko | Opera onal op miza on in district hea ng systems with the use of thermal energy storage[END_REF].

! . n ‡ˆq S . • S (X) ≤ n S (X) ≤ n ‡ˆq S . • S (X) (12)
where • S (X) is an integer variable represen ng the ON or OFF state of the equipment i and ! is the ra o between the minimum and maximum power of the equipment i. Because of the product of two decision variables, Equa on ( 12) is nonlinear and must be rewriGen using a big-M constraint formula on [START_REF] Griva | Linear and Nonlinear Op6miza6on[END_REF], introducing an intermediate variable n ‡ˆqSTX S ( ) represen ng a dynamic bound and M a very large number, as shown in Equa on [START_REF] Koch | Intelligent Controlling of Power Driven Solid Biomass CHP Plants in Flexible District Hea ng With a Seasonal Heat Storage And a Power-To-Heat Component[END_REF]. Storage Power Capacity Bounds. The power charged to or discharged from the storage at each me step is bounded by minimum and maximum values which are also decision variables, as shown in Equa on [START_REF] Dahl | Cost sensi vity of op mal sector-coupled district hea ng produc on systems[END_REF].

' 0 ≤ n N_ WX (X) ≤ n ‡ˆq WX . (1 -• WX (X)) 0 ≤ n RSWN_ WX (X) ≤ n ‡ˆq WX . • WX (X) (14)
where • WX (X) is an integer variable equal to 1 when the storage is discharging and 0 when it is charging, forbidding the storage to charge and discharge at the same me. Similarly to Equa on [START_REF] Leśko | Opera onal op miza on in district hea ng systems with the use of thermal energy storage[END_REF], Equa on ( 14) is nonlinear and must be rewriGen using a big-M constraint formula on [START_REF] Griva | Linear and Nonlinear Op6miza6on[END_REF] similar to Equa on [START_REF] Koch | Intelligent Controlling of Power Driven Solid Biomass CHP Plants in Flexible District Hea ng With a Seasonal Heat Storage And a Power-To-Heat Component[END_REF].

Storage Energy Capacity Bounds. The energy content of the storage at each me step is bounded by a maximum value which is also decision variable, as shown in Equa on [START_REF] Østergaard | Business and socioeconomic assessment of introducing heat pumps with heat storage in small-scale district hea ng systems[END_REF].

0 ≤ ^WX (X) ≤ ^ ‡ˆq WX ( 15
)
Minimum ON Time of the different equipment. In order to prevent the different equipment to start and stop too frequently, a minimum me of opera on aTer a startup has been added as a constraint (see Equa on ( 16) using a formula on similar to the one used by Yang et al. [START_REF] Yang | A novel projected two-binary-variable formula on for unit commitment in power systems[END_REF] for the so-called unit commitment problem in the electrical power produc on field.

ƒ … † • S (X) -• S (X -s) ≤ ‰ S (X) ≤ • S (X)
" . ‰ S (X) ≤ m • S (X + " -s)

= "• j @fg [START_REF] Kwon | Priority order in using biomass resources -Energy systems analyses of future scenarios for Denmark[END_REF] where " represents the minimum me ON of equipment i, and ‰ S and • S are integer decision variables already introduced beforehand and respec vely represen ng the star ng of the equipment and the state ON or OFF of the equipment.

Periodic Constraints

Addi onally, some variables such as the energy content of the storage ^WX and the ON/OFF state of the equipment • S are subjected to a periodic constraint sta ng that values at ini al and final me steps are iden cal.

Annual ε-constraints

More than only minimizing the combined CAPEX and OPEX (see Equa on ( 9)), other annual objec ves faced by the DHN operator are added. Instead of implemen ng complex and computa onally costly Mul -Objec ve Op miza on (MOO), these objec ves are added by the means of ε-Constraints as explained by Haymes et al. [START_REF] Haimes | On a bicriterion forma on of the problems of integrated system iden fica on and system op miza on[END_REF]. It is a simple approach to MOO but it requires a preselec on of the values of these objec ves prior to the op miza on which might lead to unfeasible problem. The paragraphs below present the three different ε-Constraints considered in the present study.

Maximum amount of Biomass available

As discussed in sec on 1, biomass must be seen as a limited resource. Thus, an annual ε-constraint is added as shown in Equa on ( 17)

m -SQ (X). l ≤ 8(9 < . < . < % bcdeb fg ( 17
)
where n -SQ (X) [kW] is the power supplied by the biomass boiler at each me step, < [kWhth/kg] is the Low Hea ng Value of the biomass (set to 3.8 kWhth/kg), 8(9

<

[kg] is the amount of biomass available for one year, < is the biomass boiler efficiency (set to 90%).

Maximum CO2 content of the produc on

In general, one of the main advantage of DHNs compared to other hea ng solu ons is to exhibit a low emission of CO2. In France, the building envelope regulatory framework is for example relaxed if the DHN produc on to which the building is connected is below 50g/kWth. The annual ε-constraint presented in Equa on [START_REF] Lund | Implementa on strategy for small CHP-plants in a compe ve market: the case of Lithuania[END_REF] thus forces the plant to operate below a predefined level of emission. ≤ , [START_REF] Lund | Implementa on strategy for small CHP-plants in a compe ve market: the case of Lithuania[END_REF] where [g/kWth] and , [g/kWth] are respec vely the yearly CO2 content of the plant produc on as defined respec vely in Equa on [START_REF] Vandermeulen | Controlling district hea ng and cooling networks to unlock flexibility: A review[END_REF]and the yearly CO2 content epsilon constraint of the plant produc on.

Minimum Renewable Energy Ra o (REnR) of the produc on

In France, DHN benefits from a reduced VAT if the REnR of their produc on is above 50%. Then, it is interes ng to add an annual ε-constraint on the REnR, as shown in Equa on [START_REF] Giraud | Op mal control of district hea ng systems using dynamic simula on and mixed integer linear programming[END_REF].

& '& ≥ & '& , (19)
where & '& [%] represents the yearly REnR of the plant heat produc on as defined in Equa on [START_REF]Étude de valorisa on du stockage thermique et du power-to-heat[END_REF], and & '& , [%] is the REnR epsilon constraint of the produc on.

Optimal sizing implementation and solving

The op miza on problem is implemented inside an in-house C++ based framework called PEGASE [START_REF] Vallée | An efficient co-simula on and control approach to tackle complex mul -domain energe c systems: concepts and applica ons of the PEGASE plakorm[END_REF] which embeds MILP formula on capabili es based on the Eigen linear algebra library [START_REF] Giraud | Presenta on, Valida on and Applica on of the DistrictHea ng Modelica Library[END_REF].

The overall sizing op miza on problem is solved using CPLEX [START_REF]IBM ILOG CPLEX Op6miza6on Studio V12[END_REF] for one year and a me-step of 2 hours. It leads to 35k con nuous decision variables, 18k integer decision variables and 120k constraints. The laGer leads to a computa onal me of 10 minutes for each op miza on on an office laptop. Numerous op miza ons were performed for different ε-constraints, as defined in sec on 4.2, and in different energe c contexts (French, Danish and German). The results are presented in next sec on.

4.4.Optimal sizing results

As discussed above, the focus is on the French energe c context for which the present heat produc on plant is meaningful. For comparisons purpose, the sizing is also performed with the REnR and CO2 content of the electrical produc on from Denmark and Germany (see Figure 3 and Figure 4).

Figure 5 presents the LCOE-REnR Pareto front obtained for the three countries with no constraint on the CO2 emissions and an available quan ty of biomass equivalent to 36% of the total yearly load of the DHN. It is first observed that the studied produc on plant is intrinsically highly renewable (i.e. high x-axis values). Second, the LCOE increases with an increase of the minimal REnR constraint for France and Germany but not for Denmark.

For the two former, the increase in LCOE highlights an increase of investment into a biomass boiler and a larger storage. For Denmark for which the electricity REnR is the highest (see Figure 4), the heat pump alone is sufficient even for these high constraints. Finally, it is worth men oning that, for this available biomass constraint of 36% of the total DHN demand, France cannot reach REnR higher than 85% with its current electrical mix.

Figure 6 presents the LCOE-CO2 content Pareto front for the three countries with no constraint on the available amount of biomass and a REnR constraint of 85%. That figure shows that the CO2 constraint has no effect on the French case due to its decarbonized nuclear-based electricity. It is also shown that the higher is the CO2 content of the electricity and the more it is required to invest in biomass and then in larger storage to cope with the CO2 constraint. Figure 7 presents the LCOE-REnR Pareto front obtained in France for two different levels of availability for biomass and no constraint on the CO2 content (which has no effect, see Figure 6). It shows how increasing the availability of biomass affects the Pareto front. Two main conclusions can be drawn from this figure: i) the LCOE increases almost linearly with the REnR constraint un l a transi on point at which it rapidly increases, and ii) increasing the available quan ty of biomass first displaces this transi on point towards higher REnR and second, allows reaching feasible solu on at higher REnR constraint (85% for 36% biomass available and 90% at 58% biomass available). This transi on point is linked to the need of dras cally increasing the capacity of the storage. The laGer is underlined by studying more specifically two different sizing S1 and S2, located in Figure 7. The main characteris cs of these two systems are listed in Table 3. It is shown that they are not that different in terms of heat pump and biomass capaci es. However, storages characteris cs are completely different with S2 highligh ng a clear inter-seasonal behavior compared to S1: i) 40 mes higher capacity, ii) 2 cycles (see Equa on ( 20)) instead of 48, and iii) a maximum discharge me at maximum power (see Equa on ( 21)) of more than 500 hours. That inter-seasonal behavior of S2 compared to the behavior of sizing S1 is illustrated in the yearly profiles of Figure 8

Sizing Verification through Operation

While the calcula ons that led to sizing S1 in Sec on 4 were performed with simple quasi-linear models, the purpose is now to evaluate the real performance of the produc on plant with a detailed non-linear simulator.

As explained in the methodology presented in Sec on 3, the same input data and boundary condi ons are applied to the simulator. The laGer is controlled by a MPC using a modified MILP formula on from the one of Sec on 4. Implementa on and solving of the produc on plant opera onal control are then presented. Finally, indicators of the more realis c performances are recalculated and compared to the ones obtained during the sizing-stage.

5.1.Operational MILP

The MILP formula on for the opera on is slightly different than the one presented in Sec on 4, which dealt with both the sizing and the opera on. Table 4 summarizes the differences between the sizing and the opera onal MILP formula ons.

Operational MILP implementation and solving

First, a detailed dynamic thermal-hydraulic model of the produc on plant has been implemented using Modelica programming language with the Standard Modelica and DistrictHea6ng [START_REF] Giraud | Presenta on, Valida on and Applica on of the DistrictHea ng Modelica Library[END_REF] validated libraries. In this model, all the non-linear phenomena linked to the conserva ons of mass, momentum and energy are modelled using 1D unitary models of the produc on equipment and storage. The conserva on equa ons are discretized using a finite volume scheme and relying on a staggered mesh for the fluid domain. The same approach is also used to model the "solid" domain where the heat conduction problem is solved.

For the produc on equipment, thermal iner a of the solid parts and heat losses to the ambient are respec vely modelled using equivalent, i.e. lumped, thermal capacitance and heat transfer coefficient. Regarding the heat pump, a half-Carnot model has been used for the COP.

For the sensible storage, stra fica on is accounted for with the usage of a 1D ver cal mesh for which axial conduc on and heat losses through the lateral walls are considered for each discre zed element.

Finally, to increase the me-step that can be used by the numerical integrator, the dynamics of the various PID controllers (e.g. supply temperature controller, power discharge controller …) are ignored by the simulator. This is done by formula ng the regulators model in a con nuous form. This simplifying assump ons is coherent with our objec ve to perform monthly to annual simula ons.

Second, this Modelica model is used inside an in-house simula on framework called PEGASE [START_REF] Vallée | An efficient co-simula on and control approach to tackle complex mul -domain energe c systems: concepts and applica ons of the PEGASE plakorm[END_REF] which, in addi on to the MILP formula on capabili es already men oned in Sec on 4.3, also embeds a master of cosimula on compa ble with the FMI 2.0 standard.

Instead of doing a single op miza on for a horizon of 1 year as in Sec on 4, we now realize an op miza on every 15 minutes for the next 24h. With the MPC strategy, at each me step, the MILP problem is formulated and solved using the GLPK open source solver [START_REF] Sotware Founda On | GNU Linear Programming Kit for GLPK Version 4.55[END_REF] and the simulator uses the outputs of the MPC strategy, sending back to the control module state returns such as the level of energy in the storage.

At a me step of 15 minutes, each opera onal op miza on problem has about 700 con nuous decision variables, 360 integer decision variables and 2400 constraints. On an office laptop, the laGer leads to a computa onal me of about 3 hours for a yearly simula on with an op miza on every 15 minutes and a receding horizon of 24h (96 me steps). In other words, each me step requires approximately half a second.

Sizing verification and realistic operational performances

Figure 12 presents the profiles of biomass, heat pump and storage power and DHN load obtained with the numerical simulator. These results are in the forms of 2D maps where the x-axis represents the hours of the days and the y-axis represents the days of the year. In this way, each horizontal line corresponds to a daily profile while each ver cal line corresponds to a yearly profile at a given hour of the day. For clarity reasons, the results focus on the January to May hea ng season only.

The two main conclusions to draw from these profiles are the following:

i)

The biomass boiler is used as base load throughout the opera on, since no sharp peaks of power are spoGed; ii)

The heat pump and the storage have a combined ac on to cope with the demand peaks. First, the heat pump is used during the nights at low electricity cost hours to charge the storage. Second, the storage is discharged during the demand peaks in the mornings when electricity costs are high thus reducing the opera onal costs. It is interes ng to no ce that the charging of the storage takes place just before the discharge so that heat losses in the storage are limited. We compared the power trajectories obtained during the sizing stage and the ones from Figure 12, i.e. the ones obtained during the opera onal control of the numerical simulator. While no back-up power (see Figure 1) was used during the sizing stage, 5.1% of the total heat demand was here supplied by the back-up gas boiler. The laGer highlights that the assump ons considered for the sizing MILP of Sec on 4 such as the neglec ng of temperature effects on , the stra fica on in the storage and the zero losses of the storage, led to an error of 5.1% on the energy mix.

Finally, we calculated in Table 5 the indicators listed in Sec on 3.4 for both the sizing and the opera onal stages. It is shown that the LCOE increases because of the storage losses and the use of the back-up. Also, the opera onal REnR is lower than what was expected during the sizing stage. The conclusion is that the results obtained during the sizing stage are too op mis c. The obtained differences highlight that a feedback between the opera onal results and the sizing results is s ll required and needs to be included in a future upgraded methodology. The later will be addressed in future studies.

Conclusion and Outlook

In the present study, op mal sizing and opera on of a DH produc on plant consis ng of a biomass heat-only generator, a PtoH equipment and a heat storage was inves gated in the French energe c context.

MILP formula on together with various ε-constraints were used to op mize the sizing. It was shown that the CO2 constraint had no influence but that when high REnR were targeted with limited availability of the biomass, investments in inter-seasonal storage were required. For comparisons purposes, the Danish and German electrical mix were used to perform the op mal sizing. In these cases, the CO2 constraint was proven to be the most decisive constraint.

One of the system sized in the French context was then tested in opera on using a numerical simulator and a MPC strategy. The laGer was shown to use at its best the associa on PtoH / storage by storing electricity in the form of heat when costs were low and restoring the heat at peak demands coincident with high electricity costs.

The methodology used here allows verifying the sizing by comparing indicators calculated with a simplified MILP formula on to indicators calculated with the simulated opera on of a non-linear numerical simulator, keeping iden cal boundary condi ons. It was here shown that a devia on of 5% of the energy mix was obtained which led to errors on the LCOE and REnR content of 2.6€/MWh, and 1.8% respec vely.

With the take-off of heat pump technology on district hea ng in France where a lot of medium-sized DHNs are already based on biomass, it is here proven that the combina on of both associated to some storage capabili es represents an aGrac ve solu on.

Addressing the topic of sizing in a more holis c approach is a future ambi on with the addi on of cogenera on plant, simultaneous hea ng and cooling heat pumps and cold storage to respond to electric, heat and cold demands. Regarding MPC op miza on, future studies will include mul -horizons op miza on in the methodology in order to be able i) to account for annual constraint, ii) properly op mize the opera on of an inter-seasonal storage and iii) minimize the observed devia ons between the sizing and opera onal stages.

Figure 1 :

 1 Figure 1: Studied combined Biomass and PtoH produc6on plant

7 M

 7 , it is set to 4.6 €/kW. Input data.

Figure 3 :Figure 4 :

 34 Figure 3: CO2 emissions of the electrical produc6on in g/kWhe for the year 2017 in France, Germany and Denmark. The equivalent Biomass CO2 emissions per kWhe is calculated accoun6ng for a factor of 3 ( of the heat pump) Figure 4: REnR of the electrical mix produc6on for the year 2017 in France, Germany and Denmark. The Biomass REnR is considered to be 1.

  1.4. ^ ‡ˆq WX [kWhth], n ‡ˆq WX [kW] and n ‡ˆq S

Figure 5 :

 5 Figure 5: Pareto Front LCOE-REnR for France, Denmark and Germany with no ε-constraint on the CO2 and an available mass of biomass equivalent to 36% of the total DHN load. HP only scenario is highlighted.

Figure 6 :

 6 Figure 6: Pareto Front LCOE-CO2 emission for France, Denmark and Germany with an ε-constraint on the REnR of 85% and an infinite mass of biomass available. Biomass only and HP only scenario are highlighted.

Figure 7 :

 7 Figure 7: Pareto Front LCOE-REnR for France with no ε-constraint on the CO2 and a variable available mass of biomass. S1 and S2 systems have been selected for comparisons

  (S1) and Figure 9 (S2) where the remaining discharge me at maximum power • = ^WX n ‡ˆq WX • is shown. Interes ngly, the weekly profiles of Figure 10 (S1) and Figure 11 (S2), shown for a week in March, are very similar. It is worth men oning that both sizing do not consider any back-up capacity in their investments.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8: Yearly behavior of the storage in system S1Figure9: Yearly behavior of the storage in system S2

Figure 12 :

 12 Figure 12: Biomass power, Heat Pump power, Storage power and DHN load for the January to May hea6ng season opera6on of the S1 system. The results are normalized by the peak load of the DHN (18MW). A posi6ve storage power means discharge while a nega6ve one means charge.

Table 1 :

 1 Inputs monthly total and DHW consump6on normalized respec6vely by the average monthly value of total and DHW consump6on, for an equivalent dwelling in Chambéry area for the averaged building stock[START_REF] Provent | Smart Grid Solaire Thermique: Rapport d'étude sur les réseaux de chaleur existants et les réseaux adaptés aux Eco-quar ers[END_REF]

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	TOT
	Monthly Consump6on of the equivalent dwelling (normalized by the average monthly value)
	2.1	1.7	1.5	1.0	0.5	0.2	0.2	0.2	0.2	1.0	1.5	2.0	12.0
	Monthly consump6on of DHW(normalized by the average monthly value)
	1.13	1.11	1.04	1.04	1.0	0.93	0.80	0.74	0.98	1.0	1.09	1.14	12

Table 2 :

 2 Parameter values for the produc6on plant equipment

	Parameters	Unit	Biomass Boiler	Heat Pump	Back up Gas Boiler
	N OPQR S N STUVWX S Y QT S P S S Z[ \ ]^T] X_ S	€/kWhth €/kW h -g/kWhth %	0.03 (Ref.[35]) 940 (Ref. [36]) 10 (Ref. [7]) 0.4 24 (Ref. [37]) 100	0.015 to 0.082 (Ref. [30]) 300 (Ref. [7]) 1 (Ref. [7]) 0.1 21 to 894 (variable and context-dependent, see sec on 3.3.2) 7 to 96 (variable and context-dependant, see sec on 3.3.2)	0.065 (Ref.[35]) 100 (Ref. [36]) 1 (Ref. [7]) 0 240 (Ref. [37]) 0
	3.3.2. Country-dependent data

Table 3 :

 3 Differences between the two systems of Figure7

	SIZING	S1	S2
	‡ -SQ [%^™ QˆR ] n ‡ˆq -SQ /n ™QˆR, ‡ˆq [-] n ‡ˆq šn /n ™QˆR, ‡ˆq [-] n ‡ˆq WX /n ™QˆR, ‡ˆq [-] ^ ‡ˆq WX /^ ‡ ˆq,›s [-] WX oe WX [-] • ‡ˆq WX [_ž	41% 0.34 0.46 0.49 1 48 20	58% (contraint) 0.44 0.40 0.74 40 2 537

Table 4 :

 4 Differences between the sizing and the opera6onal MILP

	MILP	Sizing		Opera/on
	Decision Variables	Sizing + Opera onal	Opera onal Only
	Annual Constraints	Yes		No
	Power Bounds Constraints	Equ. (12) and (14) non linear	Equ. (12) and (14) linear
	Objec/ve func/on Storage Loss	7 None	+ ; 5 + (; + 8(	; 5 + (; = 0.009ℎ hg (2%/ ¡¢)

Table 5 :

 5 Indicators Comparisons between sizing and opera6onal stage

	MILP	Sizing	Opera/on
	E (back-up) [%Eload]	0	5.1
	LCOE [€/MWhth]	31	33.6
	REnR [%]	85	83.2


Mokhtar Kirane
email: mkirane@univ-lr.fr

Ahmed Alsaedi

Bashir Ahmad

On Systems of Reaction Diffusion Equations with a Balance Law: The Sequel

 not only for classical systems but also for systems with either time fractional derivatives or space fractional derivatives accounting for anomalous diffusions.

1 Introduction.

In the early nineties, Martin proposed to find conditions on the reaction term f in the reaction diffusion system with a balance law:

) is supplemented with the usual boundary conditions and bounded given initial data:

The reaction term f is taken differentiable in both variables and satisfies the following requirement:

which ensures the non-negativity of the solutions u ≥ 0 and v ≥ 0 with the initial values u 0 (x) ≥ 0, v 0 (x) ≥ 0.

) can serve as a model for some chemical situations.

• the classical auto-catalytic reaction:

mU + nV (n + 1)V,
where U and V are chemical species whose local concentrations are

u = [U ] , v = [V ] .
According to the usual mass action law (an empirical law by Guldberg and Waage [START_REF] Guldberg | Etudes sur les affinités chimiques[END_REF]), the reaction term is:

f (u, v) = ku m v n , m, n ≥ 1, k > 0.
• the model of Selkov [START_REF] Sel'kov | Self-oscillations in glycolysis: 1. A simple kinetic model[END_REF]:

u t -a∆u = ν - k 1 uv γ 1 + Kv γ , v t -b∆v = -k 2 v + k 1 uv γ 1 + Kv γ ,
where u is the self inhibitor, generated with the coefficient of Hill [START_REF] Hill | The possible effects of the aggregation of the molecules of homogloebin on its dissociation curves[END_REF] γ ( usually large, often greater than 8), that is transformed in the activator v by the kinetic of Hill [START_REF] Hill | The possible effects of the aggregation of the molecules of homogloebin on its dissociation curves[END_REF]. Martin's model is obtained by setting ν = k 2 = K = 0.

• the model of Kumar [START_REF] Inamdar | Bounds on steady states for a nonsystemically autocatalysed reaction-diffusion system[END_REF]:

A -→ k 1 B-→ k 2 C
for which the decomposition of the concentration A depends exponentially on that of B (remember that the mass action law is an empirical one):

k 1 ∝ e β[B] -1, β > 0;
and f (u, v) = u n (e βv -1).

As said above, the main question posed by Martin is to ensure the boundedness of the components u and v of the solution of system (1)- [START_REF] Abdelmalek | A Lyapunov functional for a triangular reaction-diffusion system with nonlinearities of exponential growth[END_REF].

Observe that the component u is trivially bounded by a simple application of the maximum principle as the reaction term in the equation of u is nonpositive and the initial data is bounded. The question then is to show that v remains bounded whenever it is initially so.

It is a common belief that if the system of differential equations associated with system (1)-( 2) (the one obtained by dropping the diffusion terms) admits global solutions then solutions of the reaction-diffusion system will be global too. This is true for a single equation or quasi-monotone decreasing systems via the maximum principle but it is not true in general as the following simple example of Guedda and Kirane [START_REF] Guedda | Diffusion terms in systems of reaction diffusion equations can lead to blow up[END_REF] reveals it.

The system u = uv, t > 0, u(0

) = u 0 ∈ R, v = -2, t > 0, v(0) = v 0 ∈ R,
admits the global solution u(t) = u 0 e v 0 t-t 2 , v(t) = v 0 -2t, while the associated reaction diffusion system

u t -u xx = uv, x ∈ R, v t -v xx = -2, x ∈ R, subject to the initial data u(x, 0) = 1 √ 2π e -x2 2 , v(x, 0) = x 2 , admits the solution v(x, t) = x 2 , u(x, t) = 2 -1 4 2πcos( π 4 -2t) exp - x 2 2 tan π 4 -2t → ∞, t → π 8 .
One can find more examples in Churbanov [START_REF] Churbanov | An example of a reaction system with diffusion in which the diffusion leads to explosion[END_REF], Morgan [START_REF] Morgan | On a question of blow-up for semilinear parabolic systems[END_REF], Weinberger [START_REF] Weinberger | An example of blowup produced by equal diffusions[END_REF] and Souplet [START_REF] Souplet | A note on diffusion-induced blow-up[END_REF]; see also the review paper of Fila and Ninomiya [START_REF] Fila | Reaction-diffusion systems: blow-up of solutions that arises or vanishes under diffusion[END_REF].

2 The results

Classical diffusions

Martin's question has been answered positively in case of non-negative and bounded initial data for bounded Ω or not in the following situations:

2.1.1 the case of a nonlinearity of polynomial growth in v by:

In the case of bounded Ω:

-Masuda [START_REF] Masuda | On the global existence and asymptotic behavior of solutions of reaction-diffusion equations[END_REF], using the Lyapunov functional (a piece of art!):

L(t) = Ω (1 + δ(u + u 2 ))v n dx, δ > 0,
which allowed him to bound v n in the space L 1 (Ω); this is sufficient to obtain the boundedness of v. Moreover, he studied the large time behavior of the solution (u, v) and found it to be tending to a constant vector C = (c 1 , c 2 ), as t ↑ ∞, uniformly on Ω, where c j ≥ 0, j = 1, 2, satisfy

c 1 c 2 = 0, c 1 + c 2 = 1 |Ω| Ω (u 0 (x) + v 0 (x)) dx.
-Martin and Pierre retrieved Masuda's result by a duality argument [START_REF] Martin | Nonlinear Reaction Diffusion Systems[END_REF]; we don't comment more on the paper of Martin and Pierre as we will use their argument in the sequel.

A more general situation has been handled by Kouachi [36] by studying the following system:

u t -a∆u = f (u, v), (4)
v t -b∆v = g(u, v), (5)
subject to homogeneous Neumann boundary conditions and initial conditions

0 ≤ u(x, 0) = u 0 (x) ≤ M, 0 ≤ v(x, 0) = v 0 (x) ≤ N, (6)
even if the nonlinearities don't have a definite sign but satisfy the following dissipation condition:

λf (u, v) + g(u, v) ≤ 0, u, v ≥ 0, λ 1,
where f and g are of polynomial growth

f (u, v), g(u, v) ≤ C(1 + u + v) γ , u, v ≥ 0, γ ≥ 1.
Theorem 2.1. [START_REF] Kouachi | Existence of global solutions to reaction-diffusion systems via a Lyapunov functional[END_REF] Let Ω be bounded and a, b > 0, a = b, γ ≥ 1,

λ ≥ (a + b) 2 4ab m-1 ≥ 1,
where m ∈ N satisfies m > n(γ -1)/2. Then, for any u 0 , v 0 ≥ 0, bounded, problem (4)-( 5) admits a uniformly bounded in Ω × [0, +∞) global solution.

An example of nonlinearities satisfying the above conditions is given by:

f (u, v) = u p v q -u r v s , g(u, v) = u r v s -λu p v q ,
where p, q, r s ≥ 0, λ 1.

Proof. The proof of Kouachi's result is based on the Lyapunov functional

L(t) = Ω H m (u(x, t), v(x, t))dx, 0 < t < T max ,
where

H m (u, v) = m i=0 C i m K i 2 -i u i v m-i , C i m = m! i!(m -i)! , K ≥ a + b 2 √ ab , m ∈ N, λ ≥ K 2(m-1) .
Later on, Abdelmalek and Kouachi [START_REF] Abdelmalek | Proof of existence of global solutions for m-component reaction-diffusion systems with mixed boundary conditions via the Lyapunov functional method[END_REF] pushed further the study to handle a m × m ( m ≥ 2) system that generalizes Martin's system. Their system reads:

(u i ) t -a i ∆u i = f i (U ) , x ∈ Ω, t > 0, U = (u i ) m i=1 , i = 1, ..., m, (7)
subject to the general boundary conditions

λ i u i + (1 -λ i ) (u i ) ν = β i , x ∈ ∂Ω, t > 0, i = 1, ..., m, (8)
and the initial data

u i (0, x) = u 0 i (x) ≥ 0, x ∈ Ω, i = 1, ..., m, (9)
where a i are positive constants and 0

≤ λ i ≤ 1, 0 ≤ β i ∈ C 1 (∂Ω × R + ), i = 1, ..., m.
The reaction terms are supposed to satisfy the following conditions: (A1) Functions f i are continuously differentiable on R m + for all i = 1, ..., m, and satisfy f i (u 1 , ..., u i-1 , 0, u i+1 , ..., u m ) ≥ 0, for all u i ≥ 0; i = 1, ..., m, which allow the components of the solution to be nonnegative whenever they are so initially; (A2) Functions f i possess polynomial growth, i.e., for all i = 1, ..., m, there is an integer N ≥ 1 such that

|f i (U )| ≤ C 1 1 + m i=1 u i N for u i ≥ 0, i = 1, ..., m (10)
and

m-1 i=1 D i f i (U ) + f m (U ) ≤ C 2 1 + m i=1 u i (11)
for all u i ≥ 0, i = 1, ..., m, and some sufficiently large constants D i , i = 1, ..., m, C 1 and C 2 are positive and uniformly bounded functions, defined on R m + .

To render their result easy to read, let us set A ij = (a i + a j )/2 √ a i a j for all i, j = 1, ..., m, and let θ i , i = 1, ..., (m -1) , be positive constants such that

K l l > 0, l = 2, ..., m, (12)
where

K r l = K r-1 r-1 • K r-1 l -H r-1 l 2 , r = 3, ..., l, H r l = det 1≤i,j≤l
(a i,j ) i =l,...r+1

j =l-1,..r • k=r-2 Π k=1 (det k) 2 (r-k-2) r = 3, ..., l -1, K 2 l = a 1 a l l-1 Π k=1 θ 2(p k +1) 2 k • m-1 Π k=l θ 2(p k +2) 2 k positive value • l-1 Π k=1 θ 2 k -A 2 1l
and

H 2 l = a 1 √ a 2 a l θ 2(p 1 +1) 2 1 l-1 Π k=2 θ (p k +2) 2 +(p k +1) 2 k • m-1 Π k=l θ 2(p k +2) 2 k positive value • θ 2 1 A 2l -A 12 A 1l ,
and det 1≤i,j≤l

(a i,j ) i =l,...r+1 j =l-1,..r
denotes the determinant of the r square symmetric matrix obtained from the matrix (a i,j ) 1≤i,j≤m by removing the (r + 1) th , (r + 2) th , ..., l th rows and the r th , (r + 1) th , ..., (l -1) th columns. The elements of the matrix are:

a ij = a i + a j 2 θ p 2 1 1 ...θ p 2 (i-1) (i-1) θ (p i +1) 2 i ...θ (p(j-1)+1) 2 j-1 θ (p j +2) 2 j ...θ (p(m-1)+2) 2 (m-1) . ( 13
)
For the matrix )) be a solution of ( 7)-( 9) and that:

A =   
L(t) = Ω H pm (u 1 (t, x) , u 2 (t, x) , ..., u m (t, x)) dx, (14)
where

H pm (u 1 , ..., u m ) = pm p m-1 =0
...

p 2 p 1 =0 C p m-1 pm ...C p 1 p 2 θ p 2 1 1 ...θ p 2 (m-1) (m-1) u p 1 1 u p 2 -p 1 2 ...u pm-p m-1 m ,
with p m a positive integer and

C p i p j = p j ! p i !(p j -p i )! .
Then the functional L is uniformly bounded on the maximal interval of existence (0, T max ).

The boundedness of L is sufficient to show that the solution is global.

Concerning the application of Theorem 1.2, it particularly applies in particular to the following biochemical and chemical models.

Consider the scheme of reversible reactions

U 1 + U 2 k 1 k 2 U 3 , U 1 + U 4 k 3 k 4 U 5, U 2 + U 6 k 5 k 6 U 4 ; (15)
standard chemical kinetics lead to the following system of reaction-diffusion equations for the local concentrations u i of U i , i = 1, . . . , 6:

(u 1 ) t -d 1 ∆u 1 = -k 1 u 1 u 2 -k 3 u 1 u 4 + k 2 u 3 + k 4 u 5 , (u 2 ) t -d 2 ∆u 2 = -k 1 u 1 u 2 + k 2 u 3 -k 5 u 2 u 6 + k 6 u 4 , (u 3 ) t -d 3 ∆u 3 = k 1 u 1 u 2 -k 2 u 3 + k 5 u 2 u 6 -k 6 u 4 , (u 4 ) t -d 4 ∆u 4 = -k 3 u 1 u 4 + k 4 u 5 + k 5 u 2 u 6 -k 6 u 4 , (u 5 ) t -d 5 ∆u 5 = k 3 u 1 u 4 -k 4 u 5 -k 5 u 2 u 6 + k 6 u 4 , (u 6 ) t -d 6 ∆u 6 = -k 5 u 2 u 6 + k 6 u 4 . (16)
In the special case

k 5 = k 6 = 0, U 1 , U 2 , U 3 , U 4 , U 5 may represent, respec- tively, hemoglobin Hb, O 2 , HbO 2 , CO 2 and HbCO 2 .
2. The reversible chemical reaction:

µ 1 R 1 + µ 2 R 2 + ... + µ r R r k f kr ν 1 P 1 + ν 2 P 2 + ... + ν P , (17)
where R i and P i represent reactant and product species, respectively, and µ i , ν i are positive constants for each i.

For u i = [R i ], v i = [P i ] with k f , k r ,
the (non-negative) forward and reverse reaction rates, we have the following reaction diffusion system

(u i ) t -∇ • (d i ∇u i ) = µ i k r j=1 v ν j j -k f r j=1 u µ j j , i = 1, ...r, (v i ) t -∇ • (d r+i ∇v i ) = ν i k f r j=1 u µ j j -k r j=1 v ν j j , i = 1, ... . (18)
In both situations, the nonlinearities satisfy the hypotheses of the theorem.

3. The Brusselator (Prigogine and Col. [START_REF] Nicholis | Self-Organization in Non-equilibrium Systems[END_REF]) described by the following system

u t -a∆u = Au 2 v -(B + 1)u + 1, v t -b∆v = -Au 2 v + Bu,
is a model for the tri-molecular scheme

2U + V 3V, B + U V + D.
4. The model of Schnackenberg [START_REF] Schnackenberg | Simple chemical reaction systems with limit cycle behavior[END_REF] given by the system:

u t -a∆u = A -u + u γ v, v t -b∆v = B -u γ v + Bu, describes the stoichiometric reaction 2U + V 3V, A V U B,
a modification of the Brusselator by incorporating a further reaction step A V .

The question of Martin has also been answered positively in

2.1.2
The case of unbounded domain Ω = R N , and reaction terms of super polynomial growth type .

-The case of unbounded domain Ω = R N has been treated by:

-Collet and Xin [START_REF] Collet | Global existence and large time asymptotic bounds of L ∞ solutions of thermal diffusive combustion systems on R N[END_REF] using a localization of Masuda's Lyapunov function. and -Fitzgibbon, Morgan and Sanders [START_REF] Fitzgibbon | Global existence and boundedness for a class of inhomogeneous semilinear parabolic systems[END_REF] using estimates of the heat kernel.

case of Nonlinearity of sub-exponential growth in v

It has been tackled by:

-Haraux and Youkana [START_REF] Haraux | On a result of K. Masuda concerning reactiondiffusion equations[END_REF] using a modification of Masuda's Lyapunov functional in the form

L(t) = Ω (1 + δ(u + u 2
))e εv dx, δ, ε > 0 (ε small); they could handle Martin's system with the reaction term

f (u, v) = u n ψ(v) satisfying lim v→∞ 1 v log(1 + ψ(v)) = 0. ( 19
)
A nonlinearity reached by the result of [START_REF] Haraux | On a result of K. Masuda concerning reactiondiffusion equations[END_REF] is

ψ(v) = (e v γ -1), 0 < γ < 1.
Even though we do not know a practical situation leading to such nonlinearities, the mathematical result is interesting in itself.

-The case of nonlinearity of exponential growth in v has been studied in [START_REF] Herrero | Global existence for reaction-diffusion systems modelling ignition[END_REF] for the system

u t -a∆u = -f (u, v), t > 0, x ∈ R n , a > 0, v t -b∆v = f (u, v) t > 0, x ∈ R n , b > 0,
with positive and bounded initial data, where the reaction term and diffusion constants satisfy:

• f (0, v) ≥ 0, f (u, 0) ≥ 0; • f (u, v) ≤ Cφ(u)e σv , σ > 0,
for any constants C > 0, α > 0 and any continuous and positive function φ defined on [0, ∞) such that φ(0) = 0. Their contribution deals with the case b > a as the case a > b was investigated for any nonlinearity φ by Martin and Pierre [START_REF] Martin | Nonlinear Reaction Diffusion Systems[END_REF] using on an idea of Kanel [START_REF] Ya | The Cauchy problem for a system of semilinear parabolic equations with balance conditions[END_REF].

The exponential nonlinearities were suggested by Kumar [START_REF] Inamdar | Diffusive instability near Hopf bifurcation for exponentially autocatalyzed reaction-diffusion system[END_REF] for auto-catalytic reactions (see the motivation mentioned above).

It remains now an open problem to obtain the boundedness of v for nonlinearities of super exponential growth like f (u, v) = u n (e v 2 -1). Even though we do not know a definitely practical situation leading to such nonlinearities, the mathematical problem is of interest.

Nonlinearities without Limitation

In [START_REF] Kanel | Global solutions of reaction-diffusion systems with a balance law and nonlinearities of exponential growth[END_REF], Kanel and Kirane considered the triangular system

u t -a∆u = -f (u)g(v), t > 0, x ∈ Ω, ( 20
) v t -c∆u -d∆v = f (u)g(v), t > 0, x ∈ Ω, (21)
subject to the boundary conditions

u ν = v ν = 0 on ∂Ω (22)
and positive and bounded initial data, and proved the global boundedness of its solutions without any restriction on the reaction terms. Precisely, they obtained the following result.

Theorem 2.3. Let c > 0, a > d, (a -d)v 0 (x) ≥ cu 0 (x) and u 0 (x), v 0 (x) be bounded. Then there exists a nonnegative bounded solution (u, v) ((a-d)v ≥ cu) to system ( 21), ( 21), [START_REF] Hollis | On the Blow-up of Solutions to Some Semilinear and Quasilinear Reaction-diffusion Systems[END_REF].

Notice that the hypotheses of Theorem 2.1 ensures ensures the nonnegativity of u and v. Moreover, observe that the given problem is posed on a bounded domain.

Theorem 2.4. [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of reaction-diffusion equations with a balance law[END_REF] Assume thatg(v) is bounded in V(-∞), and c > d -a > 0. Then there exists a nonnegative bounded solution (u, v) to the problem ( 21)- [START_REF] Hollis | Global existence and boundedness in reaction-diffusion systems[END_REF].

We conclude this part by mentioning the result of Rebiai and Benachour [START_REF] Rebiai | Global classical solutions for reactiondiffusion systems with nonlinearities of exponential growth[END_REF].

Theorem 2.5. The solution of the system

u t -a∆u = -f (u, v), t > 0, x ∈ Ω ⊂ R n , v t -b∆v = g(u, v), t > 0, x ∈ Ω,
with a, b > 0, a = b, supplemented with the Robin boundary conditions

λ 1 u + (1 -λ 1 ) ∂u ∂ν = λ 2 v + (1 -λ 2 ) ∂v ∂ν = 0, t > 0, x ∈ ∂Ω,
and bounded initial data

u(x, 0) = u 0 (x) ≥ 0, v(x, 0) = v 0 (x) ≥ 0, x ∈ Ω,
under the hypotheses

• f (0, η) = g(0, η) = 0, g(ξ, 0) ≥ 0; • g(ξ, η) ≤ ψ(η)f (ξ, η), ψ ≥ 0, ψ ∈ C 1 (R + ); • There exists β ≥ 1 such that lim η→+∞ η β-1 ψ(η) = > 0; • g(ξ, η) ≤ ϕ(ξ)e αη β , α, β > 0,
is global provided that:

u 0 ∞ ≤ 8ab αβn(a -b) 2 , > 0.
Examples of nonlinearities satisfying the above hypotheses are:

f (u)g(v) = u(1 + v β-1 )e v β , f (u)g(v) = ue v β , m > 1.
Proof. The proof is based on the following Lyapunov functional and follows the arguments used by Barabanova in [START_REF] Barabanova | On the global existence of solutions of a reactiondiffusion equation with exponential nonlinearity[END_REF]:

L(t) = Ω δu + (M -u) -γ e αp(v+1) β (x, t) dx,
where the positive constants α, β, γ, δ and M are such that 2 and 2p > n.

β ≥ 1, u 0 ∞ < M < 2γ αβn , γ = 4ab (a -b)

Systems with a Triangular Matrix of Diffusion Coefficients

In [START_REF] Kirane | Global bounds and asymptotics for a system of reactiondiffusion equations[END_REF], Kirane initiated the study of existence of global solutions for the following triangular system with a balance law:

u t -a∆u = -uψ(v), t > 0, x ∈ Ω, v t -c∆u -d∆v = uψ(v), t > 0, x ∈ Ω, (23)
supplemented with the homogeneous Neumann boundary conditions

u ν = v ν = 0, t > 0, x ∈ ∂Ω (24)
and the initial conditions

u(0, x) = u 0 (x), v(0, x) = v 0 (x), x ∈ ∂Ω, ( 25
)
where Ω is a bounded domain of R N with a regular boundary ∂Ω, ν is the outward normal to ∂Ω, a, c, d are constants such that a > d > 0, c ≥ 0, c 2 ≤ 4ad.

The main result in [START_REF] Kirane | Global bounds and asymptotics for a system of reactiondiffusion equations[END_REF] is the following theorem.

Theorem 2.6. Let ψ satisfy the condition [START_REF] Hill | The possible effects of the aggregation of the molecules of homogloebin on its dissociation curves[END_REF],

u 0 (x), v 0 (x) ∈ L ∞ (Ω), (a -d)v 0 (x) ≥ cu 0 (x) ≥ 0.
Then the solution of the problem ( 25)-( 26) is nonnegative with v ≥ cu/(a -d) ≥ 0, and globally bounded in time. Moreover,

(u, v) → (0, v ∞ ) uniformly as t → ∞, with f (0, v ∞ ) = 0, where v ∞ = 1 |Ω| Ω (u 0 + v 0 )dx.
If N = 1, there exist T > 0 and C > 0 such that

v(t) -v ∞ ≤ Ce -ρ(t-T ) , if αλ 1 = h(v ∞ ), C(1 + t -T )e -ρ(t-T ) , otherwise, u(t) ≤ e -h(v∞(t-T ) , t ≥ T,
where α = min{a, d}, ρ = min{h(v ∞ ), dλ 1 }, λ 1 is the first positive eigenvalue of -∆ with a homogeneous Neumann boundary condition.

If N > 1, then there exist positive constants T and C such that

||(u(t), v(t) -v ∞ )|| ∞ ≤ C(1 + t -T ) -1 n+1 , t ≥ T.

Influence of the Boundary Conditions

Let us consider the system

T t = ∆u + QY f (T ), t > 0, x ∈ Ω, Y t = 1 L e ∆Y -Y f (T ), t > 0, x ∈ Ω,
modeling the reaction

A -→ B (combustion),
where T = T (x, t) is the temperature and Y = Y (x, t) is the fraction of mass, Q and Le are positive constants. Assume that f ≥ 0 is increasing and f (α) > 0, for some α > 0.

Let T and Y satisfy the boundary conditions

T ν = 0 and Y = h on x ∈ ∂Ω,
where h = is continuous and nonnegative.

Theorem 2.7. [START_REF] Avrin | Temperature growth and temperature bounds in special cases of combustion models[END_REF] If the solution of the given system exists for all time, then there exist t 1 > 0, A > 0, C such that

At + C ≤ T (., t) ∞ , t ≥ t 1 . (26)
This means that the solution, if it is global, can not be bounded.

Proof. [START_REF] Avrin | Temperature growth and temperature bounds in special cases of combustion models[END_REF] Let ϕ be the solution of the torsion problem

-∆ϕ = 1 on Ω, ϕ = 0 on ∂Ω.
It is easy to see that the function

Λ(t) = Ω (εY ϕ + δT )dx, ε, δ > 0 satisfies Λ (t) ≥ ε L Ω Y (-ϕ ν )dx.
The estimate ( 26) is obtained by integration.

The case of a full Matrix of Diffusion Coefficients

It seems that the study started with the paper of Kanel and Kirane [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of reaction-diffusion equations with a balance law[END_REF] for the following system:

u t -λ 1 ∆u -ε∆v = -uf (v), t > 0, x ∈ Ω, ( 27
)
v t -λ 2 ∆u -λ 3 ∆v = +uf (v), t > 0, x ∈ Ω, (28)
where the constants of diffusion λ 1 , λ 2 , λ 3 and ε are nonnegative satisfying the condition 4λ 1 λ 3 ≥ (ε+λ 2 ) 2 which reflects the parabolicity of the system, and the non-linearity f (v) is a smooth function satisfying uf (v) ≥ 0. Equations ( 27) and ( 28) are considered with the boundary conditions

u ν = v ν = 0, t > 0, x ∈ ∂Ω, (29)
and the initial data

u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), x ∈ Ω. ( 30
)
Let us set the following notations and hypotheses: N 1.

γ 1 := λ 3 -λ 1 -(λ 3 -λ 1 ) 2 + 4λ 2 ε /2ε; N 2. γ 2 := λ 3 -λ 1 + (λ 3 -λ 1 ) 2 + 4λ 2 ε /2ε; H1. w 10 (x) = w 1 (0, x) = v 0 (x) -γ 1 u 0 (x) ≥ 0; H2. w 20 (x) = w 2 (0, x) = v 0 (x) -γ 2 u 0 (x) ≥ 0; H3. λ 2 > λ 3 -λ 1 , 0 < ε 1; H4. f (v) = v n , n ≥ 1 is an odd integer. Theorem 2.8. Assume that γ 1 C p |λ 3 -λ 1 -λ 2 + ε| < (γ 1 + γ 2 ).
Then there exists a unique global solution (u(x, t), v(x, t)) of problem ( 27)- [START_REF] Kanel | Global existence and large time behavior of positive solutions to a reaction diffusion system[END_REF]. Moreover, the solution (u(x, t), v(x, t)) converges uniformly in Ω to a constant vector k = (k 1 , k 2 ) as t ↑ ∞, such that

k 1 ≥ 0, k 2 ≥ 0, k 1 f (k 2 ) = 0,
and

k 1 + k 2 = 1 |Ω| Ω (u 0 (x) + v 0 (x)) dx.
Here C p is the constant from the maximal regularity estimate [START_REF] Solonnikov | On boundary value problems fo linear parabolic systems of differential equations of general form[END_REF]:

u L p ((0,T )×Ω) + A p u L p ((0,T )×Ω) ≤ C p f L p ((0,T )×Ω)
of the heat equation

u (t) + A p u(t) = f (t), a.e t ∈ (0, T ), u(0) = 0, f ∈ L p ((0, T ) × Ω).
The result of [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of reaction-diffusion equations with a balance law[END_REF] attracted a significant attention and it was generalized by several authors, see, for example, [START_REF] Kouachi | Existence of global solutions to reaction-diffusion systems via a Lyapunov functional[END_REF] and [START_REF] Kouachi | Explicit Invariant Regions and Global Existence of Solutions for Reaction Diffusion Systems With a Full Matrix of Diffusion Coefficients and Nonhomogeneous Boundary Conditions[END_REF].

Let us consider the system

u t -a 11 ∆u -a 12 ∆v = f (u, v), t > 0, x ∈ Ω, ( 31
) v t -a 21 ∆u -a 22 ∆v = g(u, v), t > 0, x ∈ Ω, (32)
equipped with the boundary conditions:

λu + (1 -λ)u ν = β 1 , λv + (1 -λ)v ν = β 2 , t > 0, x ∈ ∂Ω, (33)
and the initial condition

u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), x ∈ Ω, ( 34
)
where Ω is a bounded domain of class C 1 , the constants a ij , i, j = 1, 2, are assumed to be positive and satisfy the condition

(a 12 + a 21 ) 2 ≤ 4a 11 a 22 ,
which reflects the parabolicity of the system and that implies that the eigenvalues λ 1 < λ 2 of the transposed of the matrix ((a ij )) i,j=1,2 are positive.

The initial data are supposed to be in the region

Σ =      {(u 0 , v 0 ) ∈ R 2 s.t λ 1 v 0 ≤ a 21 u 0 + a 22 v 0 ≤ λ 2 v 0 }, or {(u 0 , v 0 ) ∈ R 2 s.t λ 1 u 0 ≤ a 11 u 0 + a 12 v 0 ≤ λ 2 u 0 }.
The reaction terms f and g are assumed of class C 1 , polynomially bounded on Σ, (f (r, s), g(r, s)) ∈ Σ for all (r, s) ∈ ∂Σ, that is, r,s), for all r, s such that λ 1 s = a 21 r + a 22 s and a 21 f (r, s) + a 22 g(r, s) ≤ λ 2 g(r, s), for all r, s such that λ 2 s = a 21 r + a 22 s, and positive constants C and α > a 22 -λ 1 sufficiently close to a 22 -λ 1 ,

     λ 1 g(r, s) ≤ a 21 f (r, s) + a 22 g(
a 21 f (u, v) + Cg(u, v) ≤ C 1 (a 21 u + αv + 1), for all u, v ∈ Σ,
where C 1 is a positive constant. Kouachi [START_REF] Kouachi | Invariant regions and global existence of solutions for reactionndiusion systems with a full matrix of diffusion coeffi cients and nonhomogeneous boundary conditions[END_REF] proved the following result.

Theorem 2.9. Assume that the above conditions are satisfied and that a 11 f (r, s)+ a 22 > a 12 + a 21 . Then the system ( 31)-( 34) admits a global solution.

Proof. The proof is based on the fact that the function

L(t) = Ω H p (u(x, t), v(x, t)) dx, where H p (u, v) = p i=0 C i p θ i u i v p-i , C i p are the binomial coefficient, θ i = θ (p-i) 2 , i = 0, 1, . . . , p, θ 1 is the positive constant 2 √ a 11 + a 22 -a 12 + a 21 θ = (a 11 + a 22 ).
Finally, we want to mention the nice paper of Kuiper [START_REF] Kuiper | Positively invariant regions for strongly coupled reactionndiusion systems with a balance law[END_REF] who studied non-linear systems with full matrix of diffusion coefficients and with nonlinear boundary conditions.

In the second part of the paper, we present results concerning anomalous diffusion-reaction systems with a balance law.

Very recently, it is found that reaction-diffusion systems with time/space fractional derivative(s) describe practical situations (chemical, physical, biological, etc). In the next section, we present the results for such systems with a balance law.

The Case of Anomalous Diffusions

Hereafter, we will present results concerning a system with spatial fractional diffusions and a system with temporal fractional derivatives.

Fractional Diffusions

As reaction-diffusion systems are less documented, we will dwell a while on the tools needed for their studies.

Let the nonlocal operator (-∆) γ/2 , 0 < γ < 2 be defined by

(-∆) δ/2 u(x) = F -1 (|ξ| δ F(u)(ξ))(x), u ∈ S(R N ),
where S(R N ) is the Schwartz space and F is the Fourier transform, or equivalently by the formula

(-∆) δ/2 u(x) = C N P V R N u(x) -u(y) |x -y| N +δ dy.

Diagonal Systems

We start with the reaction anomalous diffusion system (in short (A-RDS))

u t + a(-∆) α/2 u = -f (u, v), v t + b(-∆) β/2 u = +f (u, v),
where f is assumed to be locally Lipschitz and satisfying the conditions:

f (0, v) ≥ 0, f (u, 0) ≥ 0, and |f (u, v)| ≤ L(R)|v| k + M (R), u, v ≥ 0, |u| ≤ R.
It may describe a chemical reaction taking place in a fractal medium.

We will use the notations :

Q T = R N × (0, T ), ||u|| p p = R N |u| p dx, ||u|| p p,T = Q T
|u| p dxdt, and the preliminary results.

Let S δ (t) := e -t(-∆) δ be the semigroup generated by the forward problem (for short (FP)):

Ψ t + (-∆) δ Ψ = 0, x ∈ R N , t > 0, Ψ(x, 0) = Ψ 0 (x), x ∈ R N . Lemma 1. Let 1 ≤ r ≤ p ≤ ∞ and Ψ(x) ∈ L r (R N ).
Then the solution of (FP) satisfies the estimate

||(-∆) ν/2 S δ (t)Ψ(x)|| p ≤ Ct -ν 2δ -N 2δ ( 1 r -1 p ) ||Ψ|| r , (35)
for δ > 0, ν ≥ 0.

Let ϕ satisfy the backward problem (for short (BP)) (ϕ(x, 0) = ϕ 0 (x) is not known):

-ϕ t + (-∆) β/2 ϕ = ϑ, (x, t) ∈ Q T , 0 < β < 2, ϕ(x, T ) = 0, x ∈ R N .
Lemma 2.

[61] Let 1 < q < ∞ and assume that ϑ ∈ L q (Q T ). Then problem (BP) has a unique solution ϕ such that

ϕ t ∈ L q (Q T ), (-∆) β/2 ϕ ∈ L q (Q T ).
Moreover, there exists a constant C(q, T ), independent of ϕ, such that

||ϕ|| q,T + ||(-∆) β/2 ϕ|| q,T + ||ϕ 0 || q ≤ C(q, T )||ϑ|| q,T . (36)
Lemma 3.

[58] Let (X , || . ||) be a Banach space and A be a positive operator on X . Then, for 0 ≤ ρ < σ < γ, there exists a constant C ρ,σ,γ such that, for x ∈ D(A γ ) ( the domain of A γ ),

||A σ x|| ≤ C ρ,σ,γ ||A γ x|| (σ-ρ)/(γ-ρ) ||A ρ x|| (γ-σ)/(γ-ρ) . (37)
Lemma 4. [START_REF] Lopez-Mimbela | Local time and Tanaka formula for a multitype Dawson-Watanabe super-process[END_REF] Let f α i , t > 0 be the transition density of the symmetric

α i -stable process in R N , e -t|ξ| α i = R N e i(x,ξ) f α i (t, x) dx, i = 1, 2
, where ( , ) is the inner product in R N . If α 1 ≤ α 2 , then there exists a constant K ≥ 1 such that, for every t > 0 and x ∈ R N ,

f α 2 (x, t) ≤ Kf α 1 (x, t α 1 α 2 ). (38)
Moreover, if t ≥ 1, then

f α 2 (x, t) = Kt N ( 1 α 1 -1 α 2 ) f α 1 (x, t) ≤ Kt N α 1 f α 1 (x, t), x ∈ R N . (39)
Local solutions to (A-RDS)

Theorem 3.1. (local existence). Let 0 < α, β ≤ 2. Given u 0 , v 0 ∈ C 0 (R N ), there exists a mild solution (u, v) ∈ C([0, T max [; C 0 (R N ) × C 0 (R N ))
of problem (A-RDS) satisfying the integral equations

u(t) = S α (t)u 0 - t 0 S α (t -s)f (u(s), v(s)) ds, v(t) = S β (t)v 0 + t 0 S β (t -s)f (u(s), v(s)) ds.
Moreover, lim

t-→Tmax {t + ||(u, v)|| (L ∞ ((0,t)×R N ) 2 } = +∞.
Theorem 3.2. Assume that u 0 , v 0 ∈ C 0 (R N ), u 0 , v 0 ≥ 0 a.e on R N . Let f satisfy the above conditions. Then problem (A-RDS) admits a nonnegative global solution.

Proof. As u 0 , v 0 ≥ 0, we have u ≥ 0 and v ≥ 0 and the estimate ||u|| p,T ≤ ||u 0 || p,T .

• Bound of v (duality argument).

Case 1: ( 0 < α ≤ β ≤ 2).

We have

(u + v) t + (-∆) β/2 (u + v) + (-∆) α/2 u -(-∆) β/2 u = 0.
Whereupon,

Q T (u + v)(-ϕ t + (-∆) β/2 ϕ) = R N (u 0 + v 0 )ϕ 0 + Q T u((-∆) β/2 ϕ -(-∆) α/2 ϕ).
or

Q T (u + v)ϑ = R N (u 0 + v 0 )ϕ 0 + Q T u((-∆) β/2 ϕ -(-∆) α/2 ϕ).
Using Lemma 2, we get ||ϕ|| q + ||ϕ|| q,T ≤ C||ϑ|| q,T , ||(-∆) β/2 ϕ|| q,T ≤ C||ϑ|| q,T .

Making use of inequality [START_REF] Kouachi | Existence of global solutions to reaction-diffusion systems with nonhomogeneous boundary conditions via a Lyapunov functional[END_REF] with ρ = 0, σ = α/2, and γ = β/2, we obtain

||(-∆) α/2 ϕ|| q,T ≤ ||(-∆) β/2 ϕ|| α β q,T ||ϑ|| 1-α β q,T ≤ C||ϑ|| α β q,T ||ϑ|| 1-α β
q,T ≤ C||ϑ|| q,T . Also, we have the estimates

Q T |u(-∆) β/2 ϕ| ≤ ||u|| p,T + ||(-∆) β/2 ϕ|| q,T , p + q = pq, Q T |u(-∆) α/2 ϕ| ≤ ||u|| p,T + ||(-∆) α/2 ϕ|| q,T , ≤ C(||u|| p,T + ||ϑ|| q,T ), p + q = pq.
Finally, we have

Q T (u + v)ϑ ≤ C(||u 0 + v 0 || p + ||u|| p,T + ||ϑ|| q,T + ||ϕ 0 || q ).
It then follows by duality that

||u + v|| p,T ≤ C(||u 0 + v 0 || p + ||u|| p,T + ||ϑ|| q,T + ||ϕ 0 || q ).
Therefore, for all p < ∞, the L p -norm of v remain bounded on Q Tmax . From the polynomial growth assumptionn on f , it follows that f (u, v) is also in

L p (Q Tmax ) for all p ∈ (1, ∞). For p > (N + 1)/2, we deduce that v ∈ L ∞ (Q Tmax ) and ||v|| ∞,Tmax ≤ C||f (u, v)|| p,Tmax < +∞.
Consequently, T max = +∞.

Case 2: ( 0 < β ≤ α ≤ 2). This case can be handled in the same way as Case 1.

Case of no restriction on f : ( 0 < β ≤ α ≤ 2).

Theorem 2. Assume that u 0 , v 0 ∈ C 0 (R N ), u 0 , v 0 ≥ 0 a.e on R N . Let f satisfy the conditions (C1) and (C2). Then problem (A-RDS) admits a nonnegative global classical solution.

Proof. Let S δ (t) := e -t(-∆) δ be the semigroup generated by (-∆) δ on R N and e -t|ξ| α = R N e i(x,ξ) f α (t, x) dx. Then we have u

(x, t) = S α (t)u 0 (x) -F α (x, t), v(x, t) = S β (t)v 0 (x) + F β (x, t),
where

F δ (x, t) = Q T f δ (y, s)u(x -y, t -s)f (v(x -y, t -s)) dyds, δ = α, δ = β.
As u ≥ 0, we have

0 ≤ F α (x, t) ≤ S α (t)u 0 (x) ≤ C||u 0 || ∞ . 20
Using Lemma 4 for t ≥ 1, we obtain

F β (x, t) = Q T f β (y, s)u(x -y, t -s)f (v(x -y, t -s)) dyds ≤ KT N α Q T f α (y, s)u(x -y, t -s)f (v(x -y, t -s)) dyds ≤ KT N α F α (x, t) ≤ CT N α ||u 0 || ∞ .
So v is bounded for any finite t, whereupon the solution is global.

The case of triangular system

it is shown in [START_REF] Ahmad | On nonlinear nonlocal systems of reaction diffusion equations[END_REF] that there exists a global solution for the triangular system

u t + a(-∆) α 2 u = -f (u)g(v), t > 0, x ∈ R N , (40)
v t + c(-∆) α 2 u + d(-∆) α 2 v = +f (u)g(v), t > 0, x ∈ R N , (41)
subject to the initial conditions

u(x, 0) = u 0 (x) ≥ 0, v(x, 0) = v 0 (x) ≥ 0, x ∈ R N . (42)
The nonlinerarities f : [0, +∞) -→ [0, +∞) and g : [0, +∞) -→ [0, +∞) are assumed to be locally Lipschitzian and satisfying the conditions:

(C1) f (0) = 0 ≤ g(0);

(C2) there exists a positive number k ≥ 1 such that:

g(v) ≤ Cv k for all v ≥ 0. Let H s q = H s q (R N )(s ≥ 0, 1 ≤ q ≤ ∞) be defined by H s q = H s q (R N ) = {f ∈ L q (R N )|(1 -∆) s/2 f ∈ L q (R N )}. Definition 3.1. A pair of functions (u, v) on [0, T ) × R N is called a mild solution of (40)-(41) in [0, T ) if the following assumptions hold: (i) u, v ∈ C((0, T ); L 1 (R N ) ∩ L ∞ (R N )); (ii) (u, v) satisfies the system of integral equations u(t, x) = (S α (t, •) * u 0 )(x) - t 0 S α (t -s, •) * f (u(s, •))g(v(s, •)) ds,
maximum principle [START_REF] Luchko | Maximum principle for the generalized time-fractional diffusion equation[END_REF].

To obtain the nonnegativity of v, observe that

v - c a -d u t + d(-∆) α/2 v - c a -d u = a + c -d a -d f (u)g(v).
As a + c > d, an application of the maximum principle [START_REF] Luchko | Maximum principle for the generalized time-fractional diffusion equation[END_REF] yields the nonnegativity of v.

Theorem 3.5. We have u, v ∈ C((0, T );

L 1 (R N ) ∩ L ∞ (R N )), for any 0 < T < ∞.
SThen, the solution of the system ( 41)-( 42) is global.

Remark.

One is tempted to consider the more general system

u t + a(-∆) α 2 u = -f (u)g(v), t > 0, x ∈ R N , v t + c(-∆) δ 2 u + d(-∆) γ 2 v = +f (u)g(v), t > 0, x ∈ R N .
However, the first problem that will be faced is the non-negativity of v.

The case of time fractional derivatives

In [START_REF] Alsaedi | Global existence and asymptotic behavior for a time fractional reaction-diffusion system[END_REF], a fractional time derivatives modification of the system of Martin has been studied. The system reads

c D β t u -∆u = -uf (v), t > 0, x ∈ Ω, c D β t v -d∆v = +uf (v), t > 0, x ∈ Ω, (43)
supplemented with the boundary and initial conditions

u ν (x, t) = v ν (x, t)(x, t) = 0 on ∂Ω × R + , (44) u
(x, 0) = u 0 (x), v(x, 0) = v 0 (x) in Ω, ( 45
)
where Ω is a regular bounded domain in R N (N ≥ 1) with smooth boundary ∂Ω, ν denotes the normal derivative on ∂Ω, d is a positive constant and u 0 and v 0 are nonnegative functions, c D β t , denotes the Caputo fractional derivative of order β ∈ (0, 1) defined for a differentiable function as

c D β t f (t) = 1 Γ(1 -β) t 0 f (τ ) (t -τ ) β dτ.
The main result of [START_REF] Alsaedi | Global existence and asymptotic behavior for a time fractional reaction-diffusion system[END_REF] reads: Theorem 3.6. Let u 0 , v 0 ∈ D(A) be such that u 0 and v 0 are nonnegative and bounded functions. Then the system ( 43)-( 44)-( 45) admits a unique global strong solution which satisfies u ≥ 0, v ≥ 0, u(., t) ∞ ≤ u 0 ∞ E β (-γt β ), γ > 0, t > 0, v(., t) -

1 | Ω | Ω (u 0 + v 0 ) dx ∞ ≤ CE β (-γt β ), t > 0,
where C is a positive constant and E β (z) is the Mittag-Leffler function defined by

E β (z) = +∞ k=0 z k Γ(βk + 1)
, β > 0, z ∈ C.

Proof. First, as the initial data is nonnegative, the solution is nonnegative by simple use of the maximum principle.

As uf (v) ≥ 0, we immediately get

u(., t) ∞ ≤ E β (-dt β A)u 0 ∞ ≤ u 0 ∞ ; (46)
therefore, u is uniformly bounded. For the L ∞ bounds of v, we proceed as follows. The function w = du + v satisfies

c D β t w + Aw = (d -1) c D β t u, ( 47
)
where A is the laplacian with domain

D(A) = u ∈ L 2 (Ω) / u ν = 0; Au 2 L 2 (Ω) = +∞ k=1 |λ k (u, ϕ k )| 2 < +∞ ,
where 0 = λ 0 ≤ λ 1 ≤ λ 2 ≤ .. are the eigenvalues of A and {ϕ n } n≥0 denotes the orthonormal eigenfunction system corresponding to {λ n } n≥0 ( Aϕ n = λ n ϕ n ). Applying J β t to both sides of (47), we get

J β t C D β t w + J β t Aw = (d -1)J β t c D β t u. (48)
As J β t c D β t f (t) = f (t) -f (0), c D β t J β t f (t) = f (t), t > 0, and J β t A = AJ β t , it follows that c D β t J β t w + AJ β t w = (d -1)(u -u 0 ) + w 0 .

The function W (x, t) := J β t w is then the solution of c D β t W + AW = (d -1)u + u 0 + v 0 [START_REF] Morgan | Global existence for semilinear parabolic systems[END_REF] Using ( 55) and ( 56), we obtain the estimate f (v) L q p ((0,T );L ∞ (Ω))

≤ M 1 ( Mv q ∞ T p q + C q (T )) + M 2 T p q

≤ M (1 + T p q ). ( 57
)
As u is bounded, we obtain the estimate

v ∞ ≤ v 0 ∞ + C βΓ(β) t 0 (t -s) β-1 f (v) ∞ ds.
Using the Hölder inequality and (57), it holds, for β > p q , t 0 (t -s) β-1 f (v) ∞ ds ≤ t 0 (t -s) (β-1)q q-p ds q-p q f (v) L q p ((0,t);L ∞ (Ω))

≤ M t β (1 + t -p q ).
Hence, for any t > 0, we have

v ∞ ≤ M (1 + t β ). (58)
Consequently, there is a unique globally bounded solution.

Large time estimate of u.

As c D β t v(x, t) -∆v(x, t) ≥ 0, [START_REF] Wang | Abstract fractional Cauchy problems with almost sectorial operators[END_REF] we have

v(x, t) ≥ +∞ n=0 (v 0 , ϕ n )E β (-λ n t β )ϕ n (x), (60)
where (., .) is the usual scalar product in L 2 (Ω). As lim t→+∞ E β (-λ n t β ) = 0, for λ n > 0 [START_REF] Krägeloh | Two families of functions related to the fractional powers of generators of strongly continuous contraction semigroups[END_REF], we obtain

lim t→+∞ +∞ n=0 (v 0 , ϕ n )E β (-λ n t β )ϕ n (x) = (v 0 , ϕ 0 )ϕ 0 = Ω v 0 dx. (61)
In view of ( 59) and ( 60), we obtain that for t T > 0,

v(x, t) ≥ 1 |Ω| Ω v 0 dx. (62)
So, there exists a positive constant γ such that f (v) ≥ γ, for all t T.

(63)

  Suppose that the functions f i , i = 1, ..., m, possess polynomial growth and satisfy condition[START_REF] Collet | Global existence and large time asymptotic bounds of L ∞ solutions of thermal diffusive combustion systems on R N[END_REF] for some positive constants D i , i = 1, ..., m, sufficiently large. Let (u 1 (t, .) , u 2 (t, .) , ..., u m (t, .

		a 11 a 12 . . . a 1n . . . . . . . . . . . .	   ,
		a n1 a n2 . . . a nn
	we have used the notation
	det k =	a 11 a 12 . . . a 1k . . . . . . . . . . . .	.
		a k1 a n2 . . . a kk
	Abdelmaled and Kouachi [3] proved the following result.
	Theorem 2.2.

(4i) for any T 0 (0 < T 0 < T ), sup t∈(0,T 0 )

As argued in [START_REF] Kakehi | Blowup and global existence of a solution to a semilinear reaction-diffusion system with the fractional Laplacian[END_REF], one can easily prove the following result.

If T > 0 is sufficiently small, then the system ( 40)-( 41) admits a unique mild solution

Following the lines of [START_REF] Kakehi | Blowup and global existence of a solution to a semilinear reaction-diffusion system with the fractional Laplacian[END_REF], we have the existence of strong solutions. Theorem 3.4. Assume that 1 < α < 2. Let (u, v) be the mild solution of ( 40)- [START_REF] Kuiper | Positively invariant regions for strongly coupled reactionndiusion systems with a balance law[END_REF] Proof. As (u, v) exists locally, the positivity of u follows by applying the as w 0 = du 0 +v 0 , supplemented with the initial condition W (0) := J β t w(t)| t=0 = 0 as w is continuous. Now, wee define two bounded linear operators M and P by

We clearly have

Using the estimate E β (-

, for some δ > 0, and the fact that Pu and P(u 0 + v 0 ) are bounded, we get

thanks to the estimate [START_REF] Prabhakar | A singular integral equation with a generalized Mittag-Leffler function in the kernel[END_REF]).

On the other hand, as W (0) = 0, c D β t PW = D β t PW , the function PW satisfies the equation

with the forcing term θ := (d -1)Pu + P(u 0 + v 0 ) and the initial condition J 1-β t W (0) = J 1 t w(t)| t=0 = 0. As u, u 0 and v 0 are bounded, we can assert that θ ∈ L q ((0, T ); L q (Ω)) for q > 1. In view of a maximal regularity result of Bajlekova [START_REF] Bajlekova | Fractional evolution equations in Banach spaces[END_REF], there exists a positive constant C q (T ) such that

Therefore, as D β t PW = Pw = P(du + v), we obtain

To estimate Mv, we first observe that

By writing v = Mv + Pv, we have

Consequently, we have

We then clearly have u(x, t) ≤ u 0 ∞ E β (-γt β ), for all t T.

To prove the estimate v(., t) -

we begin by noticing that

where C is a positive constant. So, it follows that

Combining (66) and (68), the desired estimate follows.

Remark.

In the light of the recent results of Kim and col. [START_REF] Kim | An L q (L p )-theory for the time fractional evolution equations with variable coefficients[END_REF], it will be interesting to study the existence of globally bounded solutions of the system c D σ t u + a(-∆)

, t > 0, x ∈ R N , where c D η t is the Caputo time fractional derivative of order 0 < η < 1.