Ilias Chorfi
email: ilias.chorfi@st.com

Corinne Alonso
email: corinne.alonso@laas.fr

Romain Monthéard
email: romain.montheard@cea.fr

Thierry Sutto
email: thierry.sutto@st.com

A GaN-Based Three-Level Dual Active Half Bridge Converter With Active Cancellation of the Steady-State DC Offset Current

Keywords: Dual active bridge (DAB), Offset current, Onboard Charger, GaN, bidirectional, multilevel

The dual active bridge (DAB) is a widely used isolated DC-DC converter topology, particularly in automotive on-board charger (OBC) applications due to its high power density, inherent bidirectional power flow and reliability. Besides, the implementation of multilevel switching cells helps reduce the voltage stress in all devices, allowing the converter to be used in high voltage battery systems (600 V and higher). Due to circuit non-idealities, a steady-state DC offset current occurs at the transformer which leads to increased power losses, loss of zero-voltage switching and possibly magnetic core saturation. This paper discusses the physical implementation of an active cancellation method of the steady-state DC offset current, applied to a GaN-based three-level dual active half bridge (DAHB). The proposed method uses a duty cycle modulation in the primary half-bridge to compensate the non-ideal behavior of the devices. Experimental validation is done on a 1.5 kW three-level dual active half bridge converter based on 650 V GaN transistors.

I. INTRODUCTION

Recently, EVs have been gaining popularity due to their contribution in reducing greenhouse gas by using cleaner energies [START_REF] Khaligh | Comprehensive Topological Analysis of Conductive and Inductive Charging Solutions for Plug-In Electric Vehicles[END_REF] and their performance [START_REF] Ahmad | A Comprehensive Review of Wireless Charging Technologies for Electric Vehicles[END_REF]. However, the range of EVs remains one of the major obstacle of the wide adoption of EV technology [START_REF] Keshmiri | Current Status and Future Trends of GaN HEMTs in Electrified Transportation[END_REF]. Reducing the weight of the power electronics while increasing the conversion efficiency can help increase the range of EVs. Furthermore, upgrading 400 V battery system to 800 V battery system helps to reduce the weight of the wiring harness and charging time [START_REF] Jung | Power Up with 800-V Systems: The benefits of upgrading voltage power for battery-electric passenger vehicles[END_REF].

Wide-bandgap (WBG) devices are expected to enable the next generation of EVs battery charging systems due to their performances such as low R DS,on and high switching frequency [START_REF] Han | Adoption of wide bandgap technology in hybrid/electric vehicles-opportunities and challenges[END_REF]. Silicon-carbide (SiC) and Gallium Nitride (GaN) materials are the best candidates for WBG devices due to their technological process maturity [START_REF] Millán | A Survey of Wide Bandgap Power Semiconductor Devices[END_REF]. GaN in particular is suitable for medium power applications.

The Dual Active bridge (DAB) is a popular choice in automotive applications in general and in DC-DC stage in OBCs in particular [START_REF] Assadi | Active Saturation Mitigation in High-Density Dual-Active-Bridge DC-DC Converter for On-Board EV Charger Applications[END_REF] due to its high power density, high efficiency, ease of zero-voltage switching (ZVS), galvanic isolation, bidirectional power flow and reliability [START_REF] De | 1.2kW dual-active bridge converter using SiC power MOSFETs and planar magnetics[END_REF]. It is based on two bridges linked through a high frequency AClink transformer, allowing bidirectional power flow. Advantages and disadvantages of different DAB configurations are presented in [START_REF] Higa | Dual active bridge DC-DC converter using both full and half bridge topologies to achieve high efficiency for wide load[END_REF]. In 800 V OBCs DC-DC stage, the utilization of 650 V GaN HEMT devices, which are readily available, require the use of multilevel topologies to lower the voltage stress across the devices down to an acceptable level. Furthermore, multilevel topologies help reduce the voltage change rate dv/dt which reduce electromagnetic interference (EMI), as well as reduce the output voltage harmonic content. In addition, employing multilevel topologies in DAB converter generates additional voltage levels, which adds more control degrees of freedom [START_REF] Gao | A DAB Converter Constructed by Nine-Switch Five-Level Active-Neutral-Point-Clamped Bridges[END_REF] [START_REF] Moonem | Analysis of a Multilevel Dual Active Bridge (ML-DAB) DC-DC Converter Using Symmetric Modulation[END_REF].

However, any mismatch between the different devices parameters, such as delay in gate signals, unequal turn-on and turn-off times and dispersion of the multiple devices R DS,on cause a steady-state DC offset current in the transformer which can lead to the magnetic core saturation.

This paper reports the physical implementation of a GaNbased three-level dual active half bridge (3-L DAHB) with active cancellation of the steady-state DC offset current. In section II, the operation of the 3-L DAHB with the dual phase shift control is explained. In addition, the DC offset current phenomenon is explained and the proposed solution in the literature are discusses. In section III, the active cancellation technique is presented and verified in simulation. The physical implementation of the 3-L DAHB converter as well as the active cancellation technique are detailed in section IV, and experimental results are discussed in section V.

II. CONTROL OF A 3-L DAHB CONVERTER

A. Modulation strategy of the 3-L ANPC leg

As shown in Fig. 1, the proposed three-level DAHB converter is made up of two three-level ANPC legs, one halfbridge on each side of the transformer. Each half-bridge consists of 6 active devices, and is able to generate a three-level voltage. S 1 -S 8 are regular NPC switching devices, while D 1 -D 4 are active clamping devices. Table I shows the switching states corresponding to the three level voltage applied to the primary side ( Vin 2 , 0 and -Vin 2 ). Phase-Shifted Pulse Width Modulation (PS-PWM) is applied to the primary and secondary ANPC half-bridges in order to generate the three-level square voltage. Fig. 2 shows the gate drive signals of the different devices in the proposed converter, all devices switch at 50% duty cycle with a phase shift applied between the outer switches (S 1 -S 4 and S 5 -S 8 ) and inner switches (S 2 -S 3 and S 6 -S 7 ). The phase shift d 1 is applied to the primary half-bridge, while d 2 is applied to the secondary half-bridge. To simplify the analysis and the design of the proposed converter, d 1 and d 2 are considered equal and are referred to as d 1 in the following. PS-PWM modulation is selected for its multiple advantages such as natural voltage balancing of the different devices, even losses distribution and ease of implementation [START_REF] Li | Generalized Phase-Shift PWM for Active-Neutral-Point-Clamped Multilevel Converter[END_REF]. The natural voltage balancing is highly desirable as it eliminates the need of using sophisticated and resource-heavy active voltage balancing algorithms.

TABLE I SWITCHING STATES OF A THREE-LEVEL ANPC HALF-BRIDGE

S 1 /D 1 S 2 S 3 S 4 /D 2 Vp 1 1 0 0 V in 2 1 0 1 0 0 0 0 1 1 -V in

B. Power flow equation

Similar to the two-level DAB converter, the power flow of the three-level DAHB converter is controlled by adjusting the phase shift Φ between the two square voltages of the primary and secondary half-bridges. This phase shift creates a voltage across the series inductor, which leads to a power transfer from the leading half-bridge to the lagging half-bridge. The transferred power from the primary side to the secondary side is described by (1) [START_REF] Moonem | Analysis of a Multilevel Dual Active Bridge (ML-DAB) DC-DC Converter Using Symmetric Modulation[END_REF], in which P out is the average output power, ω is the switching angular frequency, and n is the turn ratio of the AC-link transformer. This equation is valid only in the case of d 1 <Φ< π 2 [START_REF] Moonem | Analysis of a Multilevel Dual Active Bridge (ML-DAB) DC-DC Converter Using Symmetric Modulation[END_REF]. The value and the direction (direct or reversed) of the transferred power can be controlled by choosing the appropriate combination of d 1 and Φ.

P out 3L = n • V p • V s ω • L lk 3L • (Φ - Φ 2 π - d 2 1 π ) (1)
P out 2L = n • V p • V s ω • L lk 2L • (Φ - Φ 2 π ) (2)
If

d 1 = 0 • in (1)
, the equation is reduced to the conventional two-level DAB converter equation (2) [START_REF] Alonso | An overall study of a Dual Active Bridge for bidirectional DC/DC conversion[END_REF]. Furthermore, under the same conditions (same power, frequency and input/output voltages), the ratio between the series inductance of the twolevel DAB and three-level DAHB is always less than 1 (3). Thus, the series inductance value is smaller in three-level DAHB converter compared to the two-level DAB converter, eventually leading to greater compactness.

R = L lk 3L L lk 2L = 1 - d 2 1 π • (Φ -Φ 2 π ) , d 1 < Φ < π 2 (3)
Fig. 3 shows the different theoretical waveforms of the three-level ANPC DAHB under the condition of d 1 <Φ< π 2 . This condition can be considered the generalized case [START_REF] Moonem | Analysis of a Multilevel Dual Active Bridge (ML-DAB) DC-DC Converter Using Symmetric Modulation[END_REF].

C. DC offset current

In practice, the mismatch between the different devices along with their non-ideal behavior lead to an unbalanced voltage applied to the transformer, which causes a DC component in the input current of the transformer. This DC current can cause magnetic core saturation, which translates to an increase of the conduction losses and loss of zero voltage switching (ZVS) [START_REF] Zhang | Steady-State and Transient DC Magnetic Flux Bias Suppression Methods for a Dual Active Bridge Converter[END_REF]. Several solutions have been proposed to eliminate the problem of the steady-state DC offset current in two-level DAB converters. A straightforward method consists in adding a DC blocking capacitor in series with the inductor L lk [START_REF] Qin | A Dual Active Bridge Converter With an Extended High-Efficiency Range by DC Blocking Capacitor Voltage Control[END_REF]. However, the use of bulky series capacitors impacts considerably the power density and increases power loss. Adding an air-gap to the magnetic core decreases the permeability of the core, which makes the transformer more tolerable to the saturation [START_REF] Ortiz | Flux Balancing of Isolation Transformers and Application of "The Magnetic Ear" for Closed-Loop Volt-Second Compensation[END_REF]. Nonetheless, this method does not eliminate the DC component and leads to the decrease of the magnetizing inductance L m which increases the RMS current. In [START_REF] Ortiz | Flux Balancing of Isolation Transformers and Application of "The Magnetic Ear" for Closed-Loop Volt-Second Compensation[END_REF], it is shown that by measuring the input current or the magnetic flux of the transformer, an active feedback loop can be implemented to eliminate the DC component by controlling the duty cycle of one or the two bridges of the DAB converter. In [START_REF] Dutta | A method to measure the DC bias in high frequency isolation transformer of the dual active bridge DC to DC converter and its removal using current injection and PWM switching[END_REF], the DC offset current is measured by implementing a digital averaging filter and then used to implement a feedback PI loop on the primary bridge of the DAB converter. This method requires the sampling frequency of the averaging filter to be much higher than the switching frequency, which can be difficult to implement without a very high performance controller. The steady state DC offset current is caused by both bridges of the DAB converter. Thus, in order to eliminate completely the DC offset current, [START_REF] Zhang | Steady-State and Transient DC Magnetic Flux Bias Suppression Methods for a Dual Active Bridge Converter[END_REF] proposes a measure of the DC offset current on both sides of the transformer using Tunnel-Magnetoresistor (TMR) to implement two PI feedback loops. All the aforementioned methods are applied to the classic two-level DAB topology.

Multilevel DAB topologies have the same problem. Authors proposed to eliminate the steady state DC offset current in three-level DAB converter by adding a series capacitors on each side of the transformer [START_REF] Xuan | A Three-Level Dual-Active-Bridge Converter With Blocking Capacitors for Bidirectional Electric Vehicle Charger[END_REF]. This method increases the output voltage range. Nonetheless, it requires 200 µF high current capacitors, which considerably compromises the power density of the converter. In this paper, an active method is proposed to eliminate the steady state DC offset current without using series capacitors.

III. ACTIVE OFFSET CURRENT CONTROL TECHNIQUE

As stated before, the switches of the three-level DAHB converter under dual phase shift modulation are operating at 50%. The proposed DC offset current cancellation solution consists of implementing an average current feedback loop that changes the duty cycle of the primary switches depending on the value of the DC offset current, while the secondary leg switches always operate at 50% duty cycle. This method does not require the introduction of additional bulky DC-blocking series capacitors. Fig. 4 shows the equivalent block diagram of the proposed solution. A measurement of the DC offset current is obtained by sensing the primary current and passing it through an analog low-pass filter.

A small-signal model from the duty cycle to the average input current of the transformer is derived using the method presented in [START_REF] Song | Modeling and Analysis of 2/3-Level Dual-Active-Bridge DC-DC Converters with the Five-Level Control Scheme[END_REF]. The model helps to design a controller for eliminating the DC offset current.

A simulation model of the three-level DAHB converter is developed using PLECS to validate the proposed solution in eliminating the offset current. A positive DC offset current Fig. 6 shows the waveforms of the three-level DAHB converter with the proposed controller. Since the DC offset current is positive, the controller compensates the unbalanced voltage of the primary half bridge by increasing slightly the duty cycle. The same can be done for negative DC offset current by reducing the duty cycle.

IV. PHYSICAL IMPLEMENTATION

The realized three-level DAHB converter is modular. Each ANPC leg is built as a motherboard with three daughter boards, as seen in Fig. 7. Each daughter board implements a GaN-based two-level half-bridge optimized for very low power loop stray inductance. The daughter boards are implemented using STMicroelectronics 650 V GaN HEMTs with a typical R DS,on of 75 mΩ. A series inductor is used to emulate a 1:1 transformer (without galvanic isolation). The main parameters of the realized prototype are given in Table II. The dual phase shift modulation strategy used to control the power flow of the converter is implemented on a 32bit MCU clocked at 100 MHz. The proposed DC offset cancellation controller is implemented on the same MCU. The transformer's input current is measured using a Hall-effect current sensor. Then, the sensed current is passed through a low-pass filter to extract the DC offset current, which is used by the MCU to perform the appropriate duty cycle modulation.

V. EXPERIMENTAL RESULTS

Fig. 8 shows the measured input offset current of the transformer with and without the active cancellation technique. The offset current values are extracted for V in = V out = 600 V and for different output current values up to P = 1.5 kW. The offset current of the three-level DAHB without active cancellation was measured only up to 450 W as it reaches a significant value (1.9 A) which could saturate the magnetic core and eventually damage the converter. The offset current values for the remaining output current values are extracted using the simulation model. On the other hand, the offset current is kept under 20 mA using the implemented active cancellation technique. The implemented technique works over the entire output power range.

Fig. 9 shows the experimental waveforms of the three-level DAHB converter with the active cancellation technique for V in = V out = 600 V and P out = 1.5 kW. In the example shown, the implemented cancellation technique adjusts the duty cycle of the primary half-bridge to 54.6% to compensate for switching devices non-idealities. Furthermore, this additional control technique has no noticeable effect on the converter operation is negligible. In this work, the controlled 3-L DAHB converter achieves a peak efficiency of 98.25%. This paper has presented the design and realization of a GaN-based three-level DAHB converter for high voltage automotive OBC applications. The converter works under dualphase shift modulation and features an active control method for cancelling the steady state DC offset current. The converter design and its control have been verified in simulation. Experimental validation has been done using a 1.5 kW prototype, on which the implemented active cancellation technique has shown to work over a wide range of output power. This technique eliminates the need for bulky and expensive DCblocking capacitors, which eventually leads to a better power density.

Fig. 1 .

 1 Fig. 1. The topology of the proposed bidirectional three-level ANPC dual active half bridge.

2 Fig. 2 .

 22 Fig. 2. PS-PWM gate drive signals applied to the three-level DAHB

Fig. 3 .

 3 Fig. 3. Three-level DAHB converter waveforms

Fig. 4 .Fig. 5 .

 45 Fig. 4. Proposed controller for eliminating DC offset current in three-level DAHB

Fig. 7 .

 7 Fig. 6. 3-L DAHB converter waveforms with proposed controller

5Fig. 8 .Fig. 9 .

 89 Fig. 8. Offset current vs. output current in the three-level DAHB with and without active cancellation. V in = Vout = 600 V

TABLE II PARAMETERS

 II OF THE THREE-LEVEL DAHB CONVERTER

	Parameter		Min	Value Typ.	Max
	Input voltage	V in	600 V		800 V
	Output voltage	Vout	600 V		800 V
	Rated output power	Pout			3300 W
	transformer turns ration np/ns		1:1
	Switching frequency	fsw		100 kHz
	Series inductor	L lk		47 µH


Alaia Morell
email: alaia.morell@gmail.com

Yunne-Jai Shin

Nicolas Barrier

Morgane Travers-Trolet

Ghassen Halouani

Bruno Ernande

Bioen-OSMOSE: A bioenergetic marine ecosystem model with physiological response to temperature and oxygen

Keywords: Bioenergetic, Food web, Hypoxia, Marine ecosystem model, Phenotypic plasticity, Thermal tolerance

1 -Marine ecosystem models have been used to project the impacts of climate-induced changes in temperature and oxygen on biodiversity mainly through changes in species spatial distributions and primary production. However, fish populations may also respond to climatic pressures via physiological changes, leading to modifications in their life history that could either mitigate or worsen the consequences of climate change. 2-Building on the individual-based multispecies ecosystem model OSMOSE, Bioen-OSMOSE has been developed to account for high trophic levels' physiological responses to temperature and oxygen in future climate projections. This paper presents an overview of the Bioen-OSMOSE model, mainly detailing the new developments. These consist in the implementation of a bioenergetic sub-model that mechanistically describes somatic growth, sexual maturation and reproduction as they emerge from the energy fluxes sustained by food intake under the hypotheses of a biphasic growth model and plastic maturation age and size represented by a maturation reaction norm. These fluxes depend on temperature and oxygen concentration, thus allowing plastic physiological responses to climate change.

Introduction

1.

The development of increasingly realistic marine ecosystem models (MEMs) is needed to improve understanding and knowledge about marine ecosystems, which is one of the main challenges of the UN Decade of the Oceans [START_REF] Heymans | The Ocean Decade: A True Ecosystem Modeling Challenge[END_REF]. MEMs are end-to-end models representing ecosystems from primary production to top predators, linking the species and/or functional groups via trophic interactions. These models also account for abiotic and human activity impacts on ecosystem dynamics [START_REF] Rose | End-To-End Models for the Analysis of Marine Ecosystems: Challenges, Issues, and Next Steps[END_REF][START_REF] Steenbeek | Making spatial-temporal marine ecosystem modelling better -A perspective[END_REF][START_REF] Travers | Towards end-to-end models for investigating the effects of climate and fishing in marine ecosystems[END_REF]. MEMs are still being improved through the development of sub-models that increase their reliability in supporting ecosystem-based management [START_REF] Pikitch | Ecosystem-Based Fishery Management[END_REF][START_REF] Rose | End-To-End Models for the Analysis of Marine Ecosystems: Challenges, Issues, and Next Steps[END_REF].

The rates of ocean temperature rise and deoxygenation make urgent the development of mechanistic tools to forecast realistically their impacts from the physiology of marine organisms, to the population demographic impacts and to the consequence on marine trophic webs [START_REF] Breitburg | Declining oxygen in the global ocean and coastal waters[END_REF][START_REF] Urban | Improving the forecast for biodiversity under climate change[END_REF]. Efforts to model the temperature impacts on marine biodiversity at the ecosystem level has so far focused mainly on the bottom-up effect of temperature on the ecosystem via changes in primary production [START_REF] Lefort | Spatial and body-size dependent response of marine pelagic communities to projected global climate change[END_REF][START_REF] Moullec | An End-to-End model reveals losers and winners in a warming Mediterranean Sea[END_REF] and on the distribution shift of species according to their preferred temperature [START_REF] Albouy | From projected species distribution to food-web structure under climate change[END_REF][START_REF] Fernandes | Modelling the effects of climate change on the distribution and production of marine fishes: accounting for trophic interactions in a dynamic bioclimate envelope model[END_REF][START_REF] Moullec | An End-to-End model reveals losers and winners in a warming Mediterranean Sea[END_REF][START_REF] Serpetti | Impact of ocean warming on sustainable fisheries management informs the Ecosystem Approach to Fisheries[END_REF]. Mechanistic physiological response to temperature in MEM is modeled in size spectrum models and to our knowledge is not currently incorporated into an explicit multispecies model [START_REF] Lefort | Spatial and body-size dependent response of marine pelagic communities to projected global climate change[END_REF][START_REF] Maury | An overview of APECOSM, a spatialized mass balanced "Apex Predators ECOSystem Model" to study physiologically structured tuna population dynamics in their ecosystem[END_REF]. Although oxygen concentration is considered as a main pressure on marine biodiversity [START_REF] Laffoley | Ocean deoxygenation : everyone's problem (IUCN)[END_REF], the oxygen physiological impact on marine ecosystems is still not explicitly modeled in MEMs.

The core of recent model developments linking environmental conditions and physiological response is primarily on single-species models. These frameworks mechanistically describe life history cycles and metabolic fluxes. The response of metabolic rates to temperature is used in several frameworks [START_REF] Gillooly | Effects of size and temperature on developmental time[END_REF][START_REF] Kooijman | Dynamic Energy Budget Theory for Metabolic Organisation[END_REF] which are applied to project future population dynamics and spatial distribution under climate change scenarios. The response of metabolic rates to oxygen through its impact on ingestion [START_REF] Thomas | Effects of hypoxia on metabolic functions in marine organisms: Observed patterns and modelling assumptions within the context of Dynamic Energy Budget (DEB) theory[END_REF] has been recently introduced in the Dynamic Energy Budget framework to study the impact of hypoxia on population dynamics [START_REF] Lavaud | Modeling the impact of hypoxia on the energy budget of Atlantic cod in two populations of the Gulf of Saint-Lawrence, Canada[END_REF].

The model Bioen-OSMOSE is a new framework that mechanistically describes the emergence of life history traits through an explicit description of the underlying bioenergetic fluxes and their response to food, temperature and oxygen variation in a multispecies food web model. It has been developed from the OSMOSE framework [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]www.osmose-model.org), which is an individualbased, spatially and temporally explicit, multispecies model for regional marine ecosystems.

Designed to be possibly coupled to ocean and biogeochemical models, it includes the components of the entire ecosystem, from primary production to fish populations and human fishing activity, but the core of the model describes the dynamics of fish and macroinvertebrate species. In this paper, we provide a detailed description of the principles and equations of the Bioen-OSMOSE framework, as well as parameterization guidelines (detailed in Supporting Information). An application to the North Sea ecosystem is provided as a case study example. We then confront simulation outputs from the North Sea example to observed data to assess the consistency of the new model development and explore spatial variability in fish metabolic fluxes in response to temperature and oxygen.

Method

2.

Model description 2.1.

The Bioen-OSMOSE model (Fig. 1) represents fish individual physiological responses to temperature and oxygen variations and their consequences on fish communities in marine ecosystems. Biological unit, state variables and spatial characteristics 2.1.1.

The biological unit of the model is a school (a super-individual in individual-based modeling terminology). It is formed by individuals from the same species that are biologically identical. The state variables characterizing a school at time step belong to four categories (see Table 1 for state variable definitions and their units):

-Ontogenic state of individuals described by their age , somatic mass and gonadic mass ;

-Abundance, namely the number of individuals in the school ;

-Spatial location, i.e., the grid cell where the school is located; and -Taxonomic identity, i.e., the species to which the school belongs. Below we describe the bioenergetic sub-model that we developed to describe individual life-history and its responses to environmental variations. The individuals described with this level of detail belong to high trophic level (HTL) species, mainly fish and macroinvertebrate species.

Individual life history description 2.1.2.

Individual life history emerges from underlying bioenergetic fluxes which are described according to a biphasic growth model (Fig. 1) [START_REF] Andersen | Fish Ecology, Evolution, and Exploitation: A New Theoretical Synthesis[END_REF][START_REF] Boukal | Life-history implications of the allometric scaling of growth[END_REF][START_REF] Quince | Biphasic growth in fish I: Theoretical foundations[END_REF]. The body mass-dependent energy fluxes are allocated according to physiological tradeoffs between competing processes: maintenance, somatic growth and gonadic growth. The sexual maturation of individuals relies on the concept of maturation reaction norms that depicts how the process of maturation responds plastically to variation in body growth [START_REF] Heino | Measuring probabilistic reaction norms for age and size at maturation[END_REF][START_REF] Stearns | The Evolution of Phenotypic Plasticity in Life-History Traits: Predictions of Reaction Norms for Age and Size at Maturity[END_REF]. This combination of processes mechanistically describes how somatic growth, sexual maturation and reproduction emerge from energy fluxes sustained by food intake resulting from opportunistic sizebased predator-prey interactions.

On top of the biphasic growth model, individuals' energy mobilization and maintenance energetic costs depend on dissolved oxygen concentration and temperature so that the resulting metabolic rate (the net energy available for new tissue production) and thus somatic and gonadic growth vary with these abiotic parameters in a way that conforms to the oxygen-and capacity-limited thermal tolerance theory (OCLTT; [START_REF] Pörtner | Climate change and temperature-dependent biogeography: oxygen limitation of thermal tolerance in animals[END_REF] and more generally to thermal performance curves (TPC; [START_REF] Angilletta | Thermal Adaptation: A Theoretical and Empirical Synthesis[END_REF].

In the following description, energetic fluxes are expressed in somatic mass unit equivalents under the assumption that the ratio of energy density between somatic and gonadic tissues is independent of size. All the parameters of the bioenergetic and life-history sub-model are speciesspecific parameters except one parameter is constant across species, namely the Boltzmann constant (Table 2).

Ingestion, assimilation and mobilization 2.1.3.

For an individual in school , the ingested food at time step is described by a Holling's type 1 functional response [START_REF] Holling | The Components of Predation as Revealed by a Study of Small-Mammal Predation of the European Pine Sawfly1[END_REF] that depends on its somatic mass [START_REF] Christensen | Ecopath with Ecosim: methods, capabilities and limitations[END_REF][START_REF] Holt | Climate warming causes life-history evolution in a model for Atlantic cod (Gadus morhua)[END_REF][START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF] in two ways. First, it determines the prey biomass available to an individual of school . All other fish schools and LTL organisms (from the forcing biogeochemical model) that are present in the same grid cell are potential prey if their body size is compatible with a minimum [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]) and a maximum predator to prey size ratio based on individual total length [START_REF] Travers | Twoway coupling versus one-way forcing of plankton and fish models to predict ecosystem changes in the Benguela[END_REF] so that:

∑ { ( ) ⋂ ( )} (1)
where and are the allometric length-somatic mass relationship coefficient and exponent, respectively, is the accessibility coefficient of potential prey school to school that is essentially determined by the position in the water column of species relative to species according to their life stage, and is the biomass of prey school at time step t. The maximum possible food ingestion rate scales with the mass with a scaling exponent . The ingested food can then be written as:

( ) (2)
with the maximum ingestion rate per mass unit at exponent (or mass-specific maximum ingestion rate) of individuals in school and a multiplicative factor that depends on their life stage such that:

{ if otherwise (3)
where is the age at the end of an early-life fast-growth period (e.g., larval period or the larval and post-larval period, defined according to data availability, see Supporting Information S2) and a multiplicative factor accounting for higher mass-specific ingestion rate at this stage. A portion of the ingested food is assimilated, being lost due to excretion and feces egestion.

Reserves are not modeled in Bioen-OSMOSE: the assimilated energy is directly mobilized. The difference between assimilated and mobilized energy depends on oxygen and temperature conditions (Fig. 2). Mobilized energy , referred to as active metabolic rate in the ecophysiology literature, fuels all metabolic processes such as maintenance, digestion, foraging, somatic growth, gonadic growth, etc… The mobilization of energy relies on the use of oxygen to transform the energy held in the chemical bonds of nutrients into a usable form, namely ATP [START_REF] Clarke | Energy Flow in Growth and Production[END_REF]. In consequence, the maximum possible energy mobilized depends (i) directly on dissolved oxygen saturation that sets up an upper limit to mobilization at a given temperature and (ii) as temperature increases, on the capacity of individuals to sustain oxygen uptake and delivery for ATP production.

The mobilized energy rate is thus described by:

[ ] (4)
with [ ] and being the mobilization responses to dissolved oxygen saturation [ ] [ ] and temperature ( ), respectively, encountered by school in the grid cell . These are scaled between 0 and 1 such that, in optimal oxygen saturation and temperature conditions, all assimilated energy can be mobilized and, in suboptimal conditions, only a fraction of assimilated energy can be mobilized .

More precisely, the effect of dissolved oxygen is described by a dose-response function [START_REF] Thomas | Effects of hypoxia on metabolic functions in marine organisms: Observed patterns and modelling assumptions within the context of Dynamic Energy Budget (DEB) theory[END_REF] which increases with the saturation of dissolved oxygen: ) with parameters the asymptote and the slope of the dose-response function. The effect of temperature is such that first, energy mobilization increases with temperature according to an Arrhenius-like law due to chemical reaction rate acceleration until reaching limitation in individuals' ventilation and circulation capacity. Hence, oxygen uptake and delivery for energy mobilization saturates or even decreases at high temperatures, potentially due to temperature dependence of the rate of enzyme-catalyzed chemical reactions [START_REF] Arcus | On the Temperature Dependence of Enzyme-Catalyzed Rates[END_REF] or enzyme denaturation [START_REF] Pawar | From Metabolic Constraints on Individuals to the Dynamics of Ecosystems[END_REF]. This effect is described according to the [START_REF] Johnson | The growth rate of E. coli in relation to temperature, quinine and coenzyme[END_REF] model [START_REF] Pawar | From Metabolic Constraints on Individuals to the Dynamics of Ecosystems[END_REF]:

[ ] [ ] [ ] (5
(6)
with the Boltzmann constant, the activation energy for the Arrhenius-like increase in mobilized energy with temperature before reaching its peak value at , the activation energy when the energy mobilization declines with above , and ( ) a standardizing constant ensuring that ( ) .

Maintenance 2.1.4.

The mobilized energy fuels all metabolic processes starting in priority with the costs of maintenance of existing tissues which is often referred to as the standard metabolic rate in the ecophysiology literature. Here, we also include in the maintenance costs, the routine activities of individuals, including foraging and digestion, so that they are actually best compared to the routine metabolic rate in the ecophysiology literature. The maintenance costs are explicitly modeled to describe the share of mobilized energy between maintenance and the production of new tissues [START_REF] Charnov | Reproductive constraints and the evolution of life histories with indeterminate growth[END_REF][START_REF] Holt | Climate warming causes life-history evolution in a model for Atlantic cod (Gadus morhua)[END_REF], with precedence of the former over the latter, as well as to link mechanistically starvation mortality to energetic starvation when neither mobilized energy nor gonad energy reserves can cover the costs of maintenance (see next section on new tissue production for more details). The maintenance energy rate scales with the individual's somatic mass with the same exponent as the maximum ingestion rate. The maintenance rate also increases with the temperature experienced by individuals according to the Arrhenius law [START_REF] Brown | TOWARD A METABOLIC THEORY OF ECOLOGY[END_REF][START_REF] Gillooly | Effects of size and temperature on developmental time[END_REF][START_REF] Kooijman | Dynamic Energy Budget Theory for Metabolic Organisation[END_REF] and can be described as:

(7)

with the mass-specific maintenance rate and the Arrhenius function defined as:

(8) with the activation energy for the increase of the maintenance rate with temperature.

Net energy available for new tissue production 2.1.5.

The net energy available for new tissues production

is the difference between the mobilized energy and the maintenance costs defined as follows:

. ( 9
)
Given that the mobilized energy rate increases at a lower rate than the maintenance rate close to the species preferred temperature, it results that, all other things being equal, the emerging relationship between the net energy rate (and thus somatic and gonadic growth, see next section)

and temperature is dome-shaped and conforms to the OCLTT theory and the principle of TPC (red curve in Fig. 2). New tissue production: somatic and gonadic growth 2.1.6.

The net energy contributes to the production of new tissues with a proportion being allocated to the gonadic compartment and a proportion ( ) to the somatic one . This proportion depends on the sexual maturity status of the schools' individuals and their somatic mass . Before sexual maturation, i.e., when the maturity status , is equal to 0 and, after maturation, i.e., when , is determined such that the annual mean gonado-somatic index of individuals is constant throughout their adult life-stage and equal to [START_REF] Boukal | Life-history implications of the allometric scaling of growth[END_REF][START_REF] Lester | Interpreting the von Bertalanffy model of somatic growth in fishes: the cost of reproduction[END_REF][START_REF] Quince | Biphasic growth in fish I: Theoretical foundations[END_REF]:

̅̅̅̅ . ( 10
)
where, is the ratio of energy density between somatic and gonadic tissues, ̅̅̅ ∑ is the average net energy available per time step to individuals of school since their birth, with being the duration of a time step. Eq. 10 differs from a deterministic continuous time version of the same model [START_REF] Boukal | Life-history implications of the allometric scaling of growth[END_REF][START_REF] Lester | Interpreting the von Bertalanffy model of somatic growth in fishes: the cost of reproduction[END_REF][START_REF] Quince | Biphasic growth in fish I: Theoretical foundations[END_REF] where the current net energy would be used instead of the average According to the definition of , all net energy is allocated to somatic growth before maturation and it is shared between somatic and gonadic growth after, with the proportion allocated to gonads increasing with somatic mass [START_REF] Boukal | Life-history implications of the allometric scaling of growth[END_REF], which limits somatic growth as individuals become bigger. In case the mobilized energy cannot cover the maintenance costs , i.e., when , new tissue production is not possible and the gonadic compartment is resorbed to provide energy for sustaining maintenance. Somatic growth is then defined as follows:

{ ( ) (11)
and gonadic growth as:

{ (12)
where the second and third conditional formulas account for maintenance coverage by energy reserves contained in gonads. In the former case, gonads' energy can fully cover maintenance costs but in the latter it cannot, so that individuals undergo energetic starvation and incur additional starvation mortality (see Supporting Information S3).

Maturation 2.1.7.

Age and size at maturation vary strongly between individuals due to phenotypic plasticity. This plasticity in maturation is modeled by a deterministic linear maturation reaction norm (LMRN) that represents all the age-length combinations at which an individual can become mature [START_REF] Stearns | The Evolution of Life Histories[END_REF][START_REF] Stearns | The Evolution of Phenotypic Plasticity in Life-History Traits: Predictions of Reaction Norms for Age and Size at Maturity[END_REF]. In this framework, individuals become sexually mature when their growth trajectory in terms of body length intersects the LMRN. The maturity status of individuals of school at time step is thus described as:

{ (13)
with and the intercept and slope of the LMRN, respectively.

Reproduction 2.1.8.

Mature individuals spawn during the breeding season, then a gonad portion is used to release eggs, what is represented by a gonad portion released above 0. The sex-ratio is assumed to be 1:1 for all species and the number of eggs produced by school at time is defined as follows:

( ) (14)
with giving the time of the year at time step for a time step size of , and the mass of an egg.

At each time step of the breeding season (i.e., for which ( ) ) , new schools are produced by species , with the number of eggs, and thus individuals, per new school calculated as follows:

∑ (15)
with ∑ the total number of eggs produced by schools of species at time step , age of offspring set to 0, , their somatic mass to the mass of an egg, ,

and their gonadic mass to 0, . The new schools are released randomly depending on the specific larvae habitat map.

Mortality 2.1.9.

At each time step, a school experiences several mortality sources. The total mortality of a school is the sum of predation mortality caused by other schools, starvation mortality, fishing mortality, and additional mortalities (i.e. larval, senescence, diseases, and non-explicitly modeled predators). For a school , the starvation mortality results from the encountered food, the environmental abiotic variables and its maintenance rate. If the mobilized energy covers the maintenance costs , there is no starvation. However, if the mobilized energy is lower than the maintenance costs, the school has an energetic deficit. In this case, the gonad is used as a reserve. In case the gonad content does not cover the exceeding maintenance costs, the school faces starvation mortality proportionally to the remaining energetic deficit. The equations and details about all mortality processes are in Supporting Information S3. S3, and the data sources, references, and/or methodology to estimate these parameters are presented in the Supporting Information S6. The species spatial distributions are described by presence/absence maps, and informed per life stages (egg-larvae, juvenile, and adult) whenever information was available (Supporting Information S7). As individuals are represented in a 2D horizontal environment, a predator-prey accessibility matrix , used to determine the accessibility coefficient (see Eq. 1

and Table 2), is defined according to the vertical distribution overlap between potential predator and prey species possibly per life stage (Supporting Information S4, Table S5). The gonad portion released is estimated from the seasonality of eggs' release (see Supporting Information S2.5). The seasonality of eggs' release data are taken from the literature and presented in Supporting Information S8. The fishing mortality rates are size-dependent due to fisheries size-selectivity. The calibrated maximum fishing mortality rates are in Supporting Information S4, Table S3. The species-specific fishing selectivity curves are in Supporting Information S11. The larval and additional mortality are in Supporting Information S4, Table S3.

Forcing variables: low trophic levels and physical variables 2.2.2.

The Bioen-OSMOSE model is forced by temperature and oxygen variables and by LTL biomass fields.

The forcing data come from the regional biogeochemical model POLCOMS-ERSEM applied to the North Sea ecosystem [START_REF] Butenschön | ERSEM 15.06: a generic model for marine biogeochemistry and the ecosystem dynamics of the lower trophic levels[END_REF]. The modeled period is 2010-2019. There are five pelagic (micro-phytoplankton, diatoms, heterotrophic flagellates, micro-zooplankton, mesozooplankton) and three benthic (suspension feeders, deposit feeders and meiobenthos) LTL groups (Supporting Information S5): the biomass of the former is available in three dimensions and therefore integrated vertically, while the biomass of the latter is available in two-dimensions. Two other groups of large and very large benthos are set as homogeneous prey fields in space and time due to the absence of data and model output for these LTL groups. For the temperature and oxygen variables, their values are integrated over the 43 vertical layers of POLCOMS-ERSEM to force pelagic and demersal HTL species. Only the values in the deepest layer are used for benthic species. Monthly maps for each LTL group and temperature and oxygen variables are shown in Supporting Information S9.

Calibration 2.2.3.

The model is calibrated, i.e., parameters for which an independent estimator is unavailable are estimated, using maximum likelihood estimation based on an optimization method adapted to highdimensional parameter space, namely an evolutionary algorithm available in the package calibraR in R [START_REF] Oliveros-Ramos | calibrar: an R package for fitting complex ecological models[END_REF]. The algorithm explores the space of unknown parameters (referred to as "calibrated" in Supporting Information S4, Table S3) so as to maximize the likelihood obtained by comparing model outputs to observed data. Data used to calibrate Bioen-OSMOSE-NS are fisheries landings (ICES, 2019a), size-at-age from scientific surveys (NS-IBTS-Q1, North Sea International Bottom Trawl Survey (2010-2019), available online at http://datras.ices.dk) and estimated biomasses for assessed species [START_REF] Ices | Report of the Benchmark Workshop on Sandeel (WKSand 2016)[END_REF](ICES, , 2018a[START_REF] Ices | Report of the Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE) (3-10 May 2018[END_REF](ICES, , 2018c(ICES, , 2019b)). The discard rate of assessed species is low except for dab and plaice: the data used as landings and biomass for these species includes estimated discards from stock assessments. The biomasses estimated for stocks entirely located within the study area are directly used (herring, sandeel, sprat, sole, and whiting).

For stock with a wider distribution than the study area, the biomass data is taken proportional to April 20, 2022 20 total stock biomass according to the ratio between landings in the study area and total landings (mackerel, norway pout, plaice, saithe, cod, haddock, dab, hake). There is no biomass target value for unassessed species (horse mackerel, grey gurnard, hake, shrimp). The calibration is performed for an average state of the ecosystem for the period 2010-2019 by using observed data averaged over the period as target values (Supporting Information S10). The calibration is run using four phases with a new set of parameters to be estimated added at each phase for better convergence of the optimization: the first phase calibrates the LTL group accessibility coefficients only (Supporting Information S5), the larval mortalities are added for the second phase (Supporting Information S4, Table S3), the maximum ingestion rates are added on phase three (Supporting Information S4, Table S3), and the maximum fishing mortality rates and the additional mortality rates are added in the last phase (Supporting Information S4, Table S3).

The calibrated configuration is run for 80 years. The first 70 years is the spin-up period, a period during which the system stabilizes. The results presented hereafter are the years after the spin-up period. 28 replicates of the model are run with the same parameterization to account for Bioen-OSMOSE stochasticity.

Results and discussion

3.

In this paper, we present the Bioen-OSMOSE framework with its first application to the North Sea ecosystem, involving the coupling of the POLCOMS-ERSEM model for the physical and LTLs model with the HTLs Bioen-OSMOSE model. The North Sea trophic network has been intensively studied and modeled , either considering the whole ecosystem [START_REF] Blanchard | Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies size spectrum model[END_REF][START_REF] Cormon | Emergence of a new predator in the North Sea: evaluation of potential trophic impacts focused on hake, saithe, and Norway pout[END_REF][START_REF] Heath | Ecosystem limits to food web fluxes and fisheries yields in the North Sea simulated with an end-to-end food web model[END_REF][START_REF] Lewy | A stochastic age-length-structured multispecies model applied to North Sea stocks[END_REF][START_REF] Mackinson | An ecosystem model of the North Sea to support an ecosystem approach to fisheries management: description and parameterisation[END_REF] or part of it (including the English Channel) [START_REF] Girardin | Identification of the main processes underlying ecosystem functioning in the Eastern English Channel, with a focus on flatfish species, as revealed through the application of the Atlantis end-to-end model[END_REF][START_REF] Stäbler | Combining efforts to make maximum sustainable yields and good environmental status match in a food-web model of the southern North Sea[END_REF][START_REF] Travers-Trolet | Emergence of negative trophic level-size relationships from a size-based, individual-based multispecies fish model[END_REF]. This is the first time that the Bioen-OSMOSE model is used, i.e the OSMOSE model [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]) augmented with a mechanistic description of the emergence of life history from underlying bioenergetics and its response to temperature and oxygen seasonal and spatial variations. It is also, to our knowledge, the first application of a marine ecosystem model considering the impacts of physiologically-induced trait changes in response to food, oxygen, and temperature at the individual, population, and community levels.

Model evaluation

3.1.

The calibration procedure allowed us to estimate unknown parameters to obtain a model configuration that fairly accurately represents the North Sea ecosystem. Particular attention was paid Horizontal bars represent the first, second and third quartiles. The whiskers' extremities represent 1.5 times the interquartile space (the distance between the first and third quartile). The shrimp group was not represented in this graph, as available observed data are not sufficiently taxonomically resolved to be relevant for this functional group.

The simulated mean sizes-at-age correctly reproduce the observed ones (Fig. 4), supporting the credibility of the growth process described by the new bioenergetic sub-model. The Von-Bertalanffylike shape and the indefinite growth are two realistic properties reproduced with our model. As observed in the data, it can be noted that growth is faster during the first years of life and the sizes at older ages slowly tend to an infinite size. The simulated and observed sizes-at-age interquartile ranges overlap for almost all age classes of the species. The simulated size hierarchy between species is consistent with the observed one, which is a key expected property for a size-based model.

The simulated sizes-at-age 1 have generally the poorest fit to observed data. Size-at-age class 1 partly inherits uncertainties linked to size at hatching and to the growth rate at very early stages, notably the larval one, that is imperfectly accounted for by the multiplicative factor of maximum massspecific ingestion rate for the larvae at larval stage . In addition, growth during the first year is mainly driven by food limitation implying that the size-at-age class 1 is the result of a complex model adjustment between growth rate, competition and prey accessibility.

The variance in size-at-age differs between observed and simulated data, mainly for demersal species. The observed variance in sizes-at-age is the result of macro-environmental variations, i.e, in the abiotic environment [START_REF] Brown | TOWARD A METABOLIC THEORY OF ECOLOGY[END_REF][START_REF] Gislason | Size, growth, temperature and the natural mortality of marine fish[END_REF][START_REF] Thomas | Effects of hypoxia on metabolic functions in marine organisms: Observed patterns and modelling assumptions within the context of Dynamic Energy Budget (DEB) theory[END_REF] and food availability [START_REF] Brosset | Linking small pelagic dietary shifts with ecosystem changes in the Gulf of Lions[END_REF], micro-environmental variations, i.e., in undetectable or unaccounted for environmental conditions, and genetic variability in energy allocation inducing variability in size-at-age [START_REF] Enberg | Fishing-induced evolution of growth: concepts, mechanisms and the empirical evidence: Fishing-induced evolution of growth[END_REF]. In contrast, the variance of simulated size-at-age only results from macro-environmental variations. Thus, the species with observed variance higher than the simulated variance is because genetic and micro-environmental variances are not modeled here.

Comparison of observed and simulated age and size maturity ogives demonstrate the ability of Bioen-OSMOSE to correctly reproduce maturation patterns (Fig. 5). The simulated mean age at first maturation perfectly matches that observed for three species (cod, norway pout, and whiting). In observed data, age is given with yearly resolution. Therefore, we consider that a correct pattern is obtained for seven additional species for which the difference between simulated and observed age at maturity is less than one year (grey gurnard, haddock, hake, herring, mackerel, plaice and sole).

The worst deviation is obtained for sprat and saithe with simulated maturation occuring at later ages and larger sizes than that observed. Saithe and sprat have lower mean sizes at early ages than observed ones (Fig. 4) which can explain the late simulated maturation. However, both species have simulated maturation ages that stand within observed ranges, being lower than the upper bounds of observed maturation ages, namely 9 years for saithe in the North Sea [START_REF] Cohen | Gadiform fishes of the world (Order Gadiformes). An annotated and illustrated catalogue of cods, hakes, grenadiers and other gadiform fishes known to date[END_REF] and 4 years for sprat [START_REF] Ojaveer | Sprat, Sprattus sprattus balticus (Schn.)[END_REF] in the Baltic Sea (only mean values were reported for the North Sea population).

The use of LMRNs improves the description of the variability in the maturation process compared to the use of fixed age or size at maturity, as is most commonly done in marine ecosystem models, such as previously in the OSMOSE framework [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF] or in other models [START_REF] Audzijonyte | Atlantis: A spatially explicit end-to-end marine ecosystem model with dynamically integrated physics, ecology and socio-economic modules[END_REF]. Consequently, individuals from the same size or age class do not necessarily have the same maturity state (Fig. 5), which increases the realism of the life cycle description. For the majority of the species, the slope of simulated age and size maturity ogives is higher than the observed one, meaning that the observed maturation process is more variable than in simulations. As for size-atage, part of the observed variability in maturation is determined by genetic and/or microenvironmental variability [START_REF] Law | Fishing, selection, and phenotypic evolution[END_REF][START_REF] Van Wijk | Experimental harvesting of fish populations drives genetically based shifts in body size and maturation[END_REF] and is not modeled here. The population level: fisheries catches and species biomass 3.1.2.

Given the high statistical confidence in catch data, a greater weight was given to the corresponding likelihood component in the calibration process which resulted in reaching targets, i.e., simulated catches were within the range of observed values, for the majority of species (Fig. 6A). Plaice and dab are the two species for which the simulated catches were the farthest from their targets. Plaice and dab are two by-catch species that are largely discarded. These discards are estimated to be up to 40% for plaice and 90% for dab (ICES, 2018c). The catch data used as target here is reconstructed from the landings and discards estimates, which, in case of overestimation of discards, could explain the discrepancy between the target and the simulated values. The poor fit of plaice could also be partially explained because the migrations between the Eastern and Western English Channel stocks are not taken into account in Bioen-OSMOSE-NS (ICES, 2021).

The simulated biomasses are within acceptable ranges (Fig. 6B) and the resulting dominance ranking between species groups respects the ranking based on stock assessment estimations: small pelagic fish are the dominant group with herring as the main species. Demersal species have lower biomasses with saithe and haddock as the most abundant ones. Flatfish are the minority group in the system, which is dominated by plaice.

The simulated biomasses of mackerel and sandeel are underestimated compared to the stock assessment biomass estimates (Fig. 6B). Mackerel is a widely distributed stock in the North East Atlantic area and our biomass estimate in the North Sea (proportional to the total biomass in the North East Atlantic according to the ratio between North Sea landings and total landings) relies on the assumption of a uniform fishing effort in the assessed area. A higher fishing effort within the North Sea would lead to an overestimation of the target biomass for the area, which is credible since the North Sea is a historically heavily fished area. An ecopath model of the North Sea [START_REF] Mackinson | An ecosystem model of the North Sea to support an ecosystem approach to fisheries management: description and parameterisation[END_REF]) also estimated a lower biomass for this species (980 400 t). Underestimation for sandeel is more troublesome as it is a key forage species in the North Sea [START_REF] Engelhard | Forage fish, their fisheries, and their predators: who drives whom?[END_REF], for which the stock assessment is considered to be very detailed with 7 stocks in the area [START_REF] Ices | Report of the Benchmark Workshop on Sandeel (WKSand 2016)[END_REF].

The fact that the Bioen-OSMOSE-NS model does not describe the peculiar overwintering behavior of this species, which buries itself in sand and thus is less vulnerable to fishing and predation in winter [START_REF] Henriksen | Get up early: Revealing behavioral responses of sandeel to ocean warming using commercial catch data[END_REF], may explain the underestimation of its biomass. The sandeel may also be over-consumed by higher tropic level species in our model, indicating a missing forage species or an over-consumption of sandeel over LTL forced prey.

In Bioen-OSMOSE-NS, flatfish are represented by only the three main species of the North Sea ecosystem. However, there are other flatfish species each with low biomass levels (Scophthalmus maximus, Microstomus kitt, Scophthalmus rhombus, Platichthys flesus…) (NS-IBTS-Q1, DATRAS, [START_REF] Piet | On factors structuring the flatfish assemblage in the southern North Sea[END_REF] but whose total biomass is not negligible (NS-IBTS-Q1, DATRAS, [START_REF] Mackinson | An ecosystem model of the North Sea to support an ecosystem approach to fisheries management: description and parameterisation[END_REF]. Thus, the overestimation of plaice biomass may compensate for the absence of these other flatfish in the model that may leave an empty trophic niche.

The high biomass of the shrimp functional group, dominated in the ecosystem by the species Crangon crangon and Pandalus borealis, may seem surprising. However, as the micro-and mesozooplankton groups described by the biogeochemical model POLCOMS-ERSEM represent pelagic prey of sizes smaller than 0.5 cm only, we suggest that the shrimp functional group has a broader ecological role in the Bioen-OSMOSE-NS model by actually representing all LTL prey larger than 0.5 cm in the water column, whose biomass is critical to sustain the food web. These prey include demersal crustaceans with diel vertical migration such as Crangon crangon or Pandalus borealus as well as more pelagic species such as large amphipods (large Bathyporeia elegans) or euphausiids (Thysanoessa sp., Meganyctiphanes norvegica). For each predator species, the last size class (x-axis) includes all the larger individuals.

In Bioen-OSMOSE, the diet emerging from opportunistic predation reflects the species' relative abundances, their sizes and their spatio-temporal overlap. There are no pre-established predatorprey diet matrix in the parameterization of the model so confronting the output diets to observed ones, especially in terms of species composition, is a way to validate the model properties.

The simulated diets show patterns that are consistent with observations (Fig. 7). The model reproduces correctly observed ontogenetic diet shifts [START_REF] Timmerman | Seasonal and ontogenetic variation of whiting diet in the Eastern English Channel and the Southern North Sea[END_REF]. The prey composition shifts between pelagic early-life stages (size class 0-10 cm for fish species and 0-3 cm for the shrimp group) and the older life stages for all species. There are different emerging diet patterns depending on the predator's position in the water column. The pelagic species diet is dominated by phyto-and especially zoo-planktonic prey, which is consistent with studies on sprat and herring [START_REF] De Silva | Food and feeding habits of the herring Clupea harengus and the sprat C. sprattus in inshore waters of the west coast of Scotland[END_REF][START_REF] En Last | The food of herring, Clupea harengus, in the North Sea, 1983-1986[END_REF][START_REF] Raab | Dietary overlap between the potential competitors herring, sprat and anchovy in the North Sea[END_REF]. The benthic species diet is composed of benthic LTL groups and the shrimp functional group, similarly to results obtained by an isotopic study [START_REF] Timmerman | Plasticity of trophic interactions in fish assemblages results in temporal stability of benthic-pelagic couplings[END_REF], by plaice and sole stomach content studies for adults [START_REF] Rijnsdorp | Feeding of plaice Pleuronectes platessa L. and sole Solea solea (L.) in relation to the effects of bottom trawling[END_REF]) and for juveniles [START_REF] Amara | Feeding ecology and growth of O-group flatfish (sole, dab and plaice) on a nursery ground (Southern Bight of the North Sea)[END_REF] with smaller prey for sole than for plaice and dab of the same size [START_REF] Amara | Feeding ecology and growth of O-group flatfish (sole, dab and plaice) on a nursery ground (Southern Bight of the North Sea)[END_REF]. The demersal species have an intermediate diet composition with a high degree of piscivory for the larger fish. There is a steady increase in piscivory with size, mainly for demersal species, as shown empirically for whiting, cod, saithe and haddock in the area [START_REF] Robb | The food of five gadoid species during the pelagic O-group phase in the northern North Sea[END_REF][START_REF] Timmerman | Seasonal and ontogenetic variation of whiting diet in the Eastern English Channel and the Southern North Sea[END_REF] but not for norway pout [START_REF] Robb | The food of five gadoid species during the pelagic O-group phase in the northern North Sea[END_REF]. In addition, there is a significant part of benthic prey in the pelagic species diet, which correctly represents the strong pelagic-benthic coupling in this area [START_REF] Giraldo | Depth gradient in the resource use of a fish community from a semi-enclosed sea: Benthicpelagic coupling in fish diet[END_REF][START_REF] Timmerman | Plasticity of trophic interactions in fish assemblages results in temporal stability of benthic-pelagic couplings[END_REF]: the pelagic piscivorous fish (mackerel and horse mackerel) also feed on benthic prey which represents half of their diet [START_REF] Giraldo | Depth gradient in the resource use of a fish community from a semi-enclosed sea: Benthicpelagic coupling in fish diet[END_REF].

The physiological level: spatial pattern 3.2.

New outputs and original questions emerge from the physiological responses of metabolic rates to biotic and abiotic variables and can be explored with the Bioen-OSMOSE model. The representation of emergent spatially and seasonally varying bioenergetic fluxes is an example of the new features brought by Bioen-OSMOSE that can help improve our understanding of the relationship between temperature and ecosystem dynamics, which is crucial in the context of global warming [START_REF] Lindmark | Temperature impacts on fish physiology and resource abundance lead to faster growth but smaller fish sizes and yields under warming[END_REF]. This spatial and seasonal variability of metabolism in relation to temperature variation is often under-studied.

The simulated adult mean mass-specific net energy rate for new tissue production ̅̅̅, is the ratio between the population mean mass-specific net energy rate (see Eq. 9) and the weight at the exponent and it drives the energy allocated to growth and reproduction. It is spatially represented as an output example (Fig. 8). This spatial representation of emerging bioenergetic fluxes highlights the high variability of mean mass-specific net energy rate for three widely distributed species in the North Sea ecosystem with contrasted thermal preferences: cod, herring and dab (quoted by increasing physiological optimum temperature , defined in Fig. 2). The spatial pattern of mean mass-specific net energy is mainly explained by species thermal preferences. The species with the lowest thermal preference (cod) has a greater net energy acquisition in the northern part of the area where the water is colder on average. The opposite pattern emerges for the species with the highest April 20, 2022 29 thermal preference (dab). There is a better energy acquisition in the south where the average temperature is higher than in the north. A similar spatial pattern for growth rate was predicted as outputs of a single-species bioenergetic model for two thermophilic flatfish in the North Sea [START_REF] Teal | Bio-energetics underpins the spatial response of North Sea plaice (Pleuronectes platessa L.) and sole (Solea solea L.) to climate change[END_REF]. Herring, which has an intermediate thermal preference, exhibits a more spatially homogeneous emerging mean mass-specific net energy rate.

98 Figure 1 :

 981 Figure 1: Graphical description of the Bioen-OSMOSE model. In the Bioen-OSMOSE model trophic relationships emerge from spatio-temporal co-occurrence and size adequacy between predators and prey, the former resulting from ontogenic spatial distributions and possibly LTL biomass distribution. The life cycle emerges from the underlying bioenergetic fluxes that describe the internal processes from energy ingestion (which relies on the encountered prey) to growth, maturation and reproduction. The internal fluxes are partly driven by environmental conditions, i.e, temperature and oxygen.

Fish

  schools are distributed on a horizontal spatial grid that is composed of regular cells and that covers the geographical range of the ecosystem represented. A cell is characterized by its spatial coordinates, longitude and latitude , and several other variables: (i) the verticallydistributed values (z vertical layers) of k physico-chemical factors (such as temperature or the level of oxygen saturation (%) [ ] ) and (ii) the biomass of all low trophic level (LTL) groups (indexed by ) that are not explicitly modeled in Bioen-OSMOSE but provided as input from coupled hydrodynamic and biogeochemical models.

Figure 2 :

 2 Figure 2: Thermal responses of the bioenergetic fluxes from ingestion to tissue growth in Bioen-OSMOSE. The net energy rate dome-shaped curve (in red) conforms to the OCLTT theory and the principle of TPC. Food shortage impacts ingested energy and downstream fluxes. Hypoxia impacts mobilized energy and downstream fluxes. The maximum of the net energy rate (red curve) is called hereafter.

  ̅̅̅̅ . The averaging in a stochastic discrete time individual-based model such as Bioen-OSMOSE ensures a smooth increase of the proportion as individuals grow by dampening strong variations in and thus in due to the stochasticity of prey encounters and hence of the ingested energy .

Figure 3 :

 3 Figure 3: Case study area map. The 309 Bioen-OSMOSE-NS area is delimited 310 with the red line. The area is divided in 311 632 regular cells of 0.25° x 0.5° 312 delimited with black lines. 313

  results for indicators at different biological levels, the final result being a compromise between indicators used as target (size-at-age, catch and biomass) and emerging variables (maturation and diet).The individual level: size structure and maturity ogives 3.1.1.

Figure 4 :

 4 Figure 4: Boxplot of size-at-age per species for observed (pink) and simulated (blue) individual data.

Figure 5 :

 5 Figure 5: Age (A) and size (B) maturity ogives per species for observed (pink) and simulated (blue) data. Results are shown for species for which there is enough data to estimate and plot the observed age and size maturity ogives. Age data have a yearly resolution and size data a 2-centimeter resolution. The simulated (blue) and observed (pink) mean age at maturity are represented by vertical lines (A). The

Figure 6 :

 6 Figure 6: Fisheries catches (A) and biomasses (B), in thousand tons, per species for stock assessment estimates and simulated data. The boxplots represent the simulated data for 28 replicated

Figure 7 :

 7 Figure 7: Diet in percent of biomass eaten of prey species per size class of the predator species in cm.

  During the period described by our model(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), temperature is the main driver of spatial variability in bioenergetic fluxes for these three species. The spatial distribution of food has little or no impact on adult energy acquisition for these examples as we observe that simulated food ingestion frequently reaches the maximum at the adult stage (results not shown). Likewise, oxygen saturation has no impact on the emerging spatial pattern because oxygen saturation is not low enough to become a primary driver of bioenergetic fluxes (Vaquer-Sunyer & Duarte, 2008, Supporting Information S10).

Figure 8 :

 8 Figure 8: Spatial variability of the mean adult mass-specific net energy rate available for new tissue production, per model cell for cod, herring and dab. These three species are distributed over the whole modeled area although they have different optimum temperatures . Cod, herring and dab mean adult mass-specific net energy available rate for new tissue production averaged over the area are 7.4, 3.3, and 2.7 , respectively. Spatial variations of this bioenergetic flux can be driven by temperature, oxygen and food variation.

Table 1 :

 1 Variables and functions of the bioenergetics and life-history sub-models. ( :time step duration) 299

	Symbol	Description

Table 2 :

 2 Species-specific parameters of the bioenergetic and life-history sub-model. ( : time step duration). One parameter is constant across species,

	300
	301	namely the Boltzmann constant (	).
		Symbol	Description

from Size Maturity Age Length Key (SMALK) data (see Supporting Information S2.1) 2 from ecophysiological data (see Supporting Information S2.4) 3 from a combination of SMALK and thermal tolerance data (see Supporting Information S2.3) 4 from length-mass data 5 from fecundity data (see Supporting Information S2.5) Application to the North Sea ecosystem: Bioen-OSMOSE-NS

  2.2.

	2.2.1.	Area and explicit species

Acknowledgements

This work has been partially funded by the BiodivErsA and Belmont Forum project SOMBEE (BiodivScen programme, ANR contract n•ANR-18-EBI4-0003-01). Alaia Morell was supported by a PhD grant from Ifremer and Région Hauts-de-France. The POLCOMS-ERSEM projections were produced with funding from the European Union Horizon 2020 research and innovation programme under grant agreement No 678193 (CERES, Climate Change and European Aquatic Resources). The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM, http://www.ifremer.fr/pcdm) for providing DATARMOR storage, data access, computational resources, visualization and support services. Yunne-Jai Shin acknowledges funding support from the European Union's Horizon 2020 research and innovation programme under grant agreement No 869300 (FutureMARES) and the Pew marine fellows programme.

Data Archiving

The Bioen-OSMOSE-NS configuration and its associated version of the Bioen-OSMOSE model executable will be deposited on Zenodo. Model code will be available on Github. The scripts developed to estimate Bioen-OSMOSE-NS parameters are available on Github.

Applying Bioen-OSMOSE to the North Sea allows demonstrating the feasibility of its parameterization for several species with different levels of available knowledge and allows evaluating the framework capabilities. Bioen-OSMOSE-NS simulates many different outcomes that convincingly reproduce observations such as biomasses, catches, sizes at age, maturation ogives, and diets. The model also produces compelling spatial responses of the bioenergetic fluxes to temperature variations.

The Bioen-OSMOSE framework is also intended to be used for hindcast or forecast simulations.

Hindcasting could help disentangling the effects of temperature increase and/or oxygen depletion on the historical trends in life-history traits. Hindcasting with Bioen-OSMOSE could also be useful to understand the contribution of temperature-and oxygen-induced physiological changes in population and community dynamic alterations that were observed in past periods. Given the increasing need to reliably forecast biodiversity under future climate change scenarios, we believe that Bioen-OSMOSE will also allow improving projections of regional ecosystem dynamics by taking into account future individual-level physiological changes and their consequences at the population and community levels.

Conflicts of Interest

The authors declare that they do not have personal interest that could have appeared to influence the work reported in this paper.

Author Contributions

Yunne-Jai Shin and Bruno Ernande conceived and supervised the project. Bruno Ernande and Alaia