Clement Mayet
email: clement.mayet@lecnam.net

Pablo Arboleya

Alain Bouscayrol

Bassam Mohamed

Philippe Delarue

Islam El-Sayed

Different Models for the Simulation of DC Railway Systems Supplied by Non-Reversible Substations

Keywords: DC Traction System, Non-reversible Substation, Dynamic and Quasi-Static Models, Energetic Macroscopic Representation (EMR), Modified Nodal Analysis. I

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Nowadays, transportation systems are key players to face current challenges in terms of energy saving and greenhouse gas emissions [START_REF] Barrero | Energy savings in public transport[END_REF]. Public transportation systems, such as electric railways, are currently employed in urban areas due to their high transportation capacities and because they do not produce local pollutant emission [START_REF] Eastham | Running off the rails[END_REF]. However, the global warming and increasing cost of energy lead to still develop more efficient systems. Innovative solutions have been introduced to save more braking energy, such as on-board or wayside Energy Storage Systems (ESS) [START_REF] Pastor | Flywheels store to save: Improving railway efficiency with energy storage[END_REF]- [START_REF] Gee | Analysis of trackside flywheel energy storage in light rail system[END_REF] and reversible Traction Power Substations (TPS) [START_REF] Tzeng | Electric network solutions of DC transit systems with inverting substations[END_REF], [START_REF] Rufer | A supercapacitor-based energy storage substation for voltage compensation in weak transportation network[END_REF]. However, these solutions induce supplementary cost and their effectiveness must be assessed before their real implementations.

Several simulation tools have been proposed in the last decade to study railway systems and to estimate the potential energy that can be saved [START_REF] Cai | Iterative techniques for the solution of complex DC-rail-traction systems including regenerative braking[END_REF]- [START_REF] Arboleya | DC railway simulation including controllable power electronics and energy storage devices[END_REF]. Conventional TPS (nonreversible), that connect the AC grid to the DC traction network, and the squeeze control (over-voltage protection) play an important role in the solving procedure of railway simulation tools. Indeed, the conventional TPS induces non-linear behavior when it switches between on-and off-states [START_REF] Mayet | Comparison of different EMR-based models of Traction Power Substations for energetic studies of subway lines[END_REF]. The squeeze control of a braking train is activated when the supplied DC voltage is too high. In such case, the braking energy is dissipated in the rheostatic brake instead to be recovered on the DC network and induces energy discontinuities from the DC network point of view. Several simplifications on the models are commonly accepted due to the complexity of such systems.

A specific scenario has been identified as particularly difficult to study with most of the existing simulation tools due to the widely accepted simplifications. This scenario consists of a low receptive line with most of the trains in braking mode. All the TPSs are thus blocked, which disconnect the railway system from the AC grid. It becomes an islanded system in which the braking trains supply the other consumers. In such case, some simulation tools neglect what happens in the system during this period because no energy is provided by the AC grid [START_REF] Gee | Analysis of trackside flywheel energy storage in light rail system[END_REF], [START_REF] Chymera | Modeling electrified transit systems[END_REF]. Other simulation tools use complex iterative methods to approach the mathematical solution [START_REF] Cai | Iterative techniques for the solution of complex DC-rail-traction systems including regenerative braking[END_REF], [START_REF] Arboleya | DC railway simulation including controllable power electronics and energy storage devices[END_REF]. The inability of most of the simulation tools to deal with the case where all the TPSs are blocked can conduct to approximations in the design of potential innovative subsystems. Indeed, their sizing is sensitive because they are based on a cost-benefit calculation based on potential energy savings. However, these innovative subsystems are especially used when the line receptivity is low to save the braking energy, which would be dissipated in the conventional railway system. Furthermore, the control of ESS or reversible TPS is generally based on the DC network voltage [START_REF] Wang | A novel energy management strategy of onboard supercapacitor for subway applications with permanent-magnet traction system[END_REF], [START_REF] Gee | Analysis of trackside flywheel energy storage in light rail system[END_REF], [START_REF] Tzeng | Electric network solutions of DC transit systems with inverting substations[END_REF], [START_REF] Arboleya | DC railway simulation including controllable power electronics and energy storage devices[END_REF]. It is therefore crucial to be able to accurately estimate this voltage.

The objective of the paper is to compare different modelling methodologies to simulate DC railway systems supplied by non-reversible TPSs. The periods when the line receptivity is low, which conduct to block all the TPSs, are especially studied. The paper is then organized as follows. The studied case is presented in section II. Section III describes the models of the different subsystems of railway systems. Section IV is focused on different modelling approaches to simulate the entire systems. The simulation results are then compared and analyzed in Section V.

II.

DESCRIPTION OF THE STUDIED CASE

A specific scenario has been defined in order to analyze the simulation results and compare the different simulation approaches. A conventional light rail DC system with a 6.8 km line composed of 7 passenger stations (PS) and 3 conventional TPSs is considered (Fig. 1 and Fig. 2). The 3 TPSs are supposed identical and have been previously defined in [START_REF] Mayet | Comparison of different EMR-based models of Traction Power Substations for energetic studies of subway lines[END_REF]. The vehicles are light automatic subways, which have been described in [START_REF] Mayet | Comparison of different models and simulation approaches for the energetic study of a subway[END_REF]. The studied scenario considers the headway of 120s, which requires 10 trains in circulation for a simulation time of 1000s. Furthermore, it is assumed that there is no traffic initially on the line. The different trains (T1 to T10) are thus injected on the line from PS1 and PS7 (Fig. 3). The velocity profile has been defined based on the respect of the kinematical limitations (torque, speed, acceleration, voltage, etc.).

III. MODELS OF THE BASIC RAILWAY SUBSYSTEMS

This section describes the models of the basic subsystems of railway systems (i.e. TPS, Train, DC network). These models are organized using the Energetic Macroscopic Representation (EMR). The EMR is a graphical description tool, which highlights the energy properties of a system [START_REF] Bouscayrol | Systemic design methodologies for electrical energy systems -Chapter 3: Graphic formalism for the control of multi-physical energetic systems: COG and EMR[END_REF]. It organizes the system into interconnected basic elements. All elements are connected according to the interaction principle: the product of the action and reaction variables leads to the power exchanged by connected elements. Furthermore, all the components are described respecting the physical integral causality.

A. Model of the Conventional Traction Power Substations

Neglecting harmonics, the equivalent circuit (Fig. 4) of a conventional TPS is described using a switched-model to estimate its on-and off-states behaviors (Fig. 5 and Fig. 6) [START_REF] Mayet | Comparison of different EMR-based models of Traction Power Substations for energetic studies of subway lines[END_REF]. The on-state model imposes to the DC network the rectified open-circuit voltage ess0(t) through the resistance Rss (TPS Joule losses). It estimates the 3-phases AC currents ig-3φ(t) supplied by the AC grid according to the DC current iss(t) and the resistance Rir-3φ (transformer iron losses). The off-state model imposes currents equal to zero on the TPS DC (iss(t) = 0) as well as AC (it-3φ(t) = 0) sides. Switching conditions are defined (Fig. 5). Note that some models increase the value of the resistance Rss in off-state instead of using a current-voltage switching model [START_REF] Gee | Analysis of trackside flywheel energy storage in light rail system[END_REF], [START_REF] Cai | Iterative techniques for the solution of complex DC-rail-traction systems including regenerative braking[END_REF]. With such an approach, the choice of the value of Rss in offstate can impact the computation time or the simulation results themselves. More explanations about this TPS model are in [START_REF] Mayet | Comparison of different EMR-based models of Traction Power Substations for energetic studies of subway lines[END_REF].

B. Model of the Trains

The train model considers the different components of the train [START_REF] Mayet | Comparison of different models and simulation approaches for the energetic study of a subway[END_REF]. It determines the electric power Pref(t) absorbed or recovered on the train input filter by the traction and auxiliary subsystems. The auxiliary power Paux(t) is generally assumed known. The traction power Ptrac(t) is estimated according to the line profiles (slope, speed, etc.), the mass of the vehicle, the resistive efforts (friction, aerodynamic, etc.), the components losses, and the electro-kinematical limitations (torques/currents, acceleration, speed, voltage, etc.). More details can be found in [START_REF] Mayet | Comparison of different models and simulation approaches for the energetic study of a subway[END_REF]. The train equivalent circuit (Fig. 7) consists in a current source iref(t), which represents the traction and auxiliary subsystems, connected in parallel with the rheostatic brake and the DC bus capacitor Cf. The brake is activated according to the squeeze control (overvoltage protection), which induces a braking current iR(t) depending on iref(t) and uf(t) (Fig. 8). All the braking energy is therefore dissipated in the brake when the maximal voltage U2 is reached. The braking coefficient kR(t) is increased linearly from U1 to U2. Finally, the resistance Rf represents the losses in the filter as described in [START_REF] Mayet | Comparison of different models and simulation approaches for the energetic study of a subway[END_REF]. The input filter model and the squeeze control have huge impacts on the simulation approach and the solving procedure of the entire railway system. It will be discussed in Section IV.

Net. vg-3φ(t) ig-3φ(t) it-3φ(t) iir-3φ(t) Grid ess0(t) iss(t) = 0 vg-3φ(t) vg-3φ(t) vg-3φ(t) it-3φ(t) = 0 vg-3φ(t) it-3φ(t)

Model of the DC Traction Network

The model of the DC traction network assumes linear resistance distribution along the line. The resistance of each line segment Rnet(t) represents the resistance of the third rail (or catenary) plus the return rail (or ground) between two subsystems (TPS or vehicle). It is a position-dependent parameter that varies with the movement of the trains.

IV. SPECIFIC DC RAILWAY SIMULATION APPROACHES

Three specific approaches for the simulation of complete DC railway systems are described. An example is considered with 1 TPS and 2 trains to illustrate the different approaches (Fig. 9). The 3 approaches use the well-known Modified Nodal Analysis (MNA) to solve the DC traction network [START_REF] Ho | The modified nodal approach to network analysis[END_REF], [START_REF] Acary | Nonsmooth modelling and simulation of switched circuits[END_REF]. The MNA is generally expressed by [START_REF] Barrero | Energy savings in public transport[END_REF], where G is the conductance matrix of the DC network (with Gnet-ik = Rnet-ik -1

), R is the matrix that consider the different input resistances (i.e. Rss-i or Rf-i) of the voltage source subsystems, B is a matrix corresponding to the Kirchhoff currents equations and is composed of 1 for the voltage source nodes whose branch relations are introduced. The vectors J and E include respectively variables of known initial currents and voltages based on the outputs of the train and TPS models (current or voltage sources). Then, the vectors V and I are the unknown variables, respectively the nodes voltages and the currents of the voltage sources subsystems. These unknown variables are obtained by inverting equation [START_REF] Barrero | Energy savings in public transport[END_REF] and are used to define the reaction variables of the different models of the subsystems, which are simulated on the traction network.

Note that the TPS model is the same for the 3 different approaches. In on-state, it corresponds to a voltage source, which imposes the voltage ess0(t). The unknown variable iss(t) is obtained from I (MNA results). In off-state, it corresponds to a current source, which imposes a current Ji = 0. The unknown variable ess(t) is obtained from V. One can note that the MNA formulation is strongly dependent on the models of the different subsystems. Indeed, the chosen models define the inputs and outputs of the MNA, and thus the formulation. More explanations will be provided in the final paper.

Rnet-12(t) ess-1(t) iss-1(t) ess0-1(t)

Rss-1

TPS-1 N1 uf-1(t) Rf-1 it-1(t) ut-1(t) N2 TRAIN-1 Rnet-23(t) uf-2(t) Rf-2 it-2(t) ut-2(t)

N3

TRAIN-2 Fig. 9. Example of railway system with 1 TPS and 2 trains.

[ 𝐺 -𝐵 𝐵 𝑡 -𝑅 ] [ 𝑉 𝐼 ] = [ 𝐽 𝐸 ], then [ 𝑉 𝐼 ] = [ 𝐺 -𝐵 𝐵 𝑡 -𝑅 ] -1 [ 𝐽 𝐸 ] (1)

A. Dynamic model (DYN)

The first simulation approach is based on a dynamic model of the train input filter (Fig. 7) (Fig. 11.a). The DC bus imposes the voltage uf(t) as inputs for the DC network and the brake, traction, and auxiliary subsystems. The vehicle is thus a voltage source m connected to a node i. The internal DC bus current if(t) is calculated by the train model according to the currents iref(t) (Fig. 7) and iR(t) (imposed by the squeeze control in Fig. 8). Furthermore, the DC bus model imposes to the MNA the voltage uf(t). The corresponding unknown variable it(t) (current exchanged between the train and the DC network) is obtained from the MNA results with the element Im (vector I). More explanations and equations will be provided in the final paper.

B. Quasi-Static model 1 (QS1)

The QS1 model can be considered as representative of the conventional approach used in most of the existing simulation tools [START_REF] Gee | Analysis of trackside flywheel energy storage in light rail system[END_REF], [START_REF] Chymera | Modeling electrified transit systems[END_REF]. This model neglects the DC bus and considers the train as a pure current source it(t) (Fig. 11.b), with it(t) equal to if(t). Furthermore, iR(t) and iref(t) are calculated from the voltage value uf(t) of the previous simulation step-time. The current it(t) is thus imposed as a known variable in the MNA vector J. The corresponding unknown ut(t) is obtained from the MNA results (V). The actual DC bus voltages uf(t) is then estimated with [START_REF] Eastham | Running off the rails[END_REF].

𝑢 𝑓 (𝑡) = 𝑢 𝑡 (𝑡) -𝑅 𝑓 𝑖 𝑡 (𝑡) (2)
One drawback of this approach, contrary to the dynamic model-based approach, is the inability to deal with the case where all the TPSs are blocked. In such case, there is no voltage reference in the DC network, which is disconnected from the AC grid. It is therefore assumed that no energy is exchanged within the entire system: the voltages all over the DC network are imposed to the rated value and the currents are imposed to zero. More details will be provided in the final paper.

C. Quasi-Static model 2 (QS2)

The QS2 model is an intermediary approach that considers current-voltage switching sources for the train model (Fig. 11.c). The current iR(t) is null (it(t) = iref(t)) if the brake is not activated (off). In such case, the train is a current source and imposes it(t) to the MNA. The unknown variable uf(t) is obtained from the MNA results and (2). If the brake is activated (on), the train is considered as a voltage source m, which regulates the voltage to its maximal value U2 and imposes uf(t) = U2 to the MNA. The unknown variable it(t) is obtained from Im (vector I) and the current iR(t) is calculated with [START_REF] Pastor | Flywheels store to save: Improving railway efficiency with energy storage[END_REF].

𝑖 𝑅 (𝑡) = 𝑖 𝑡 (𝑡) -𝑖 𝑟𝑒𝑓 (𝑡) (3)
This approach neglects the DC bus and the squeeze control between U1 and U2 (Fig. 8). Other switching conditions are thus used to activate the brake (Fig. 10) (current/voltage sources). More details will be provided in the final paper.

V. ANALYSIS OF THE SIMULATION RESULTS

The results obtained by the 3 approaches are compared and analyzed according to the case study (section II) with the reference power Pref(t) = Ptrac(t) + Paux(t) (Fig. 12.a). The results in terms of energy in the TPSs are quite similar. The authors therefore focus in this paper on the source of errors, which is concentrated when all TPSs are blocked. In this situation, the DC traction system is completely isolated from the AC grid and there is no bus that actuates as slack. However, energy exchanges are still possible within the DC network. First, the periods when all TPSs are blocked are identified and are then referred to as "blocked periods". The powers provided by all TPSs with the dynamic model are presented (Fig. 12.b). The black signal is added to highlight these periods. All methods present a similar number of instants (between 9 and 10% of the simulation time) in which all TPSs are blocked at the same time.

The third rail or catenary voltage ut(t) and the current it(t) of Train 6 are analyzed to highlight the differences between the 3 approaches (Fig. 13). A zoom is provided for a blocked period. The DYN model (blue) gives accurate estimations even during the blocked periods. DYN has been experimentally validated in [START_REF] Mayet | Comparison of different EMR-based models of Traction Power Substations for energetic studies of subway lines[END_REF] and is considered as reference in this paper. QS1 is not able to estimate voltages and power flows during the blocked periods (10% of the time) because it returns the trivial solution (no currents on the line, and DC voltages equal to the rated value). QS2 results (black) give better results than QS1 because it can estimate the energy recovery and the voltages values during the blocked periods. However, there is still an error in comparison to the DYN simulation due to the squeeze control simplification.

An energetic study quantifies the impacts of each approach on the estimation of the internal power flows during the blocked periods. It can be noted that these power flows have no impact on the energy consumption on the AC grid because the railway system is disconnected from the grid during the blocked periods. Nevertheless, it has an impact on the DC network behavior and its different subsystems. Table I summarizes for each method the total reference energy Eref (from Pref(t)) and actual energy Ecat injected or absorbed from the catenary by all the trains. The positive (absorbed energy E (+)) and negative (recovered energy E (-)) are separated to have clear analysis on the requested energy and the regeneration capacity of all trains. The reference energy is the same for the 3 approaches. The difference between Eref and Ecat corresponds to the losses in the filters and the brakes. The errors on the absorbed and recovered energies with QS1 (vs DYN) are important during the blocked periods (100%) because the solver returns only the trivial solution. The differences are also not negligible during the entire simulation, with errors of 4% and 14.8%, respectively on the absorbed and recovered energy by all trains. The energetic error for QS2 is negligible on the entire simulation and is 4% during the blocked periods. Fig. 14 presents a specific instant (670s) where all TPSs are blocked (QS1 is not presented). The horizontal axis corresponds to the positions. The vertical axis represents the voltages at the junction with the catenary (right scalesolid red line) and the reference power Pref(t) and actual power exchanged with the DC network (left scaleblack and red bars respectively). The sum of all power references of all trains in that specific instant is -637kW. There is a power surplus that leads to increasing the voltages and blocking all TPSs. The net power injection at this specific instant is 23.4kW with DYN and 16.1kW with QS2 (0kW with QS1). Regarding the voltage profile, QS2 provides a higher voltage profile due to the squeeze control simplification, which explains the lower power injection with this method. The accuracy of the estimation of the DC voltages and the regeneration capacity of the trains is crucial to perform a study assessment of innovative subsystems. Indeed, these new subsystems are particularly used when the line receptivity is low. In addition, their controls are often based on the voltage value. It is thus essential to develop accurate and fast simulation tools. More analysis will be provided in the final paper.

VI. CONCLUSION

The paper presents and compares 3 different methodologies for the simulation of DC railway systems. The study is focused on the periods in which all the traction power substations are blocked at the same time since the results in other conditions are quite similar. A first simulation approach, based on an experimentally validated dynamic model (DYN), is presented and used as reference for comparisons. Quasi-static method 1 (QS1) provides the trivial solution and leads to significant errors compared to DYN. Quasi-static method 2 (QS2) provides better results than QS1 but needs to be further improved to obtain more accurate solutions closer to the DYN method. Future research will thus focus on improving the accuracy of the QS2 method to obtain equivalent results to those obtained with the DYN method. In addition, the potential reduction of the computation time with this new method will be studied.

Fig. 1 .

 1 Fig. 1. Electrical description of the railway system (3 TPSs).

Fig. 2 .

 2 Fig. 2. Profile description of the studied line.

Fig. 3 .

 3 Fig. 3. Schedule of the different trains.

Fig. 11 .

 11 Fig. 11. EMR of the vehicle: a) Dynamic model, b) Quasi-static model 1, c) Quasi-static model 2.

Fig. 12

 12 Fig. 12. a) Train reference power (Pref(t) = Ptrac(t) + Paux(t)), b) TPSs powers. Fig. 13. Third rail or catenary voltage and current of Train 6 with the 3 models.

Fig. 14 .

 14 Fig. 14. Simulation results at the instant 670s: a) DYN results, b) QS2 results, c) differences between DYN and QS2.

TABLE I .

 I ENERGETIC STUDY OF POWER FLOWS WITHIN THE DC TRACTION NETWORK

		REF (kWh)	DYN (kWh)	QS1 (kWh)	QS2 (kWh)	QS1 vs DYN (%)	QS2 vs DYN (%)
		Eref (+)	Eref (-)	Ecat (+)	Ecat (-)	Ecat (+)	Ecat (-)	Ecat (-)	Ecat (-)	Ecat (+)	Ecat (-)	Ecat (+)	Ecat (-)
	Entire simulation (1000s)	452.5	-142.9	469.7	-130.8	450.8	-111.4	469.6	-130.7	4.02	14.83	0.02	0.07
	All TPSs blocked periods	19.0	-28.9	19.4	-19.8	0.0	-0.0	18.6	-18.9	100.00 100.00	4.12	4.54


Mohammed Lotfi Hachemi

Abdelghani Ghomari

Yassine Hadjadj-Aoul

Gerardo Rubino

Mobile traffic forecasting using a combined FFT/LSTM strategy in SDN networks

Keywords: SDN, LSTM, time-series, forecasting, Fast Fourier Transform

Over the last few years, networks' infrastructures are experiencing a profound change initiated by Software Defined Networking (SDN) and Network Function Virtualization (NFV). In such networks, avoiding the risk of service degradation increasingly involves predicting the evolution of metrics impacting the Quality of Service (QoS), in order to implement appropriate preventive actions. Recurrent neural networks, in particular Long Short Term Memory (LSTM) networks, already demonstrated their efficiency in predicting time series, in particular in networking, thanks to their ability to memorize long sequences of data. In this paper, we propose an improvement that increases their accuracy by combining them with filters, especially the Fast Fourier Transform (FFT), in order to better extract the characteristics of the time series to be predicted. The proposed approach allows improving prediction performance significantly, while presenting an extremely low computational complexity at run-time compared to classical techniques such as Auto-Regressive Integrated Moving Average (ARIMA), which requires costly online operations.

I. Introduction

The continued emergence of new and constrained services is one of the major challenges facing Infrastructures' Providers (InPs) [START_REF] Camacho | Emerging technologies and research challenges for intelligent transportation systems: 5G, HetNets, and SDN[END_REF]. In order to cope with these new trends, InPs are currently shifting their hardware-based equipment to much more agile software-based solutions. This has been made possible through the separation of the control plane from the data plane, with the Software Defined Networking (SDN) paradigm, and the decoupling of network functions from their hardware, with the Network Function Virtualization (NFV) concept [START_REF] Wood | Toward a software-based network: integrating Software Defined Networking and Network Function Virtualization[END_REF].

Although these technologies enable InPs to meet their objectives, they nevertheless introduce new challenges [START_REF] Han | Network Function Virtualization: Challenges and opportunities for innovations[END_REF]. One of the major ones is the scaling up (respectively down) of functions to accommodate an increase (respectively decrease) in the workload. Indeed, scaling-up network functions involves the deployment of one or more new instances, which generally requires several seconds or even minutes. As a consequence and in order to avoid any congestion, several studies in the literature focus on load prediction to anticipate its possible increase [START_REF] Alawe | Improving Traffic Forecasting for 5G Core Network Scalability: A Machine Learning Approach[END_REF].

Forecasting time-series evolution is a fairly classical field of research [START_REF][END_REF] with a vast associated literature. Initial approaches were generally based on exponential smoothing or Auto-Regressive Integrated Moving Average (ARIMA) models [START_REF] Yoo | Time-series forecast modeling on high-bandwidth network measurements[END_REF]. Although effective, these techniques have recently been superseded by deep learning-based strategies [START_REF] Siami-Namini | A comparison of arima and lstm in forecasting time series[END_REF]. Among them, Recurrent Neural Networks (RNN) are certainly the most effective in recognizing traffic patterns, thanks to their recurrent structure. In particular, Long Short Term Memory (LSTM) networks outperform the majority of existing techniques and do not present the problem of the vanishing gradient from which RNNs typically suffer [START_REF] Sak | Long Short-Term Memory recurrent neural network architectures for large scale acoustic modeling[END_REF].

In this paper, we propose to go beyond existing work by putting forward a new approach to the problem of time series forecasting. The idea is to combine an LSTM network with a feature extraction module in order to facilitate the task of identifying patterns in the considered time series. In particular, we analyze the capacity of the Discrete Fourier Transform (DFT) [START_REF] Briggs | The DFT: An Owners' Manual for the Discrete Fourier Transform[END_REF] to extract the most relevant features from the time series, while canceling noise-like signals, which may affect learning. This combination not only decreases significantly the prediction error, but can also reduce the number of required neurons, thus accelerating the learning process. In addition, it maintains pretty good performance over longer prediction horizons, as can be seen in the sequel.

The rest of this paper is organized as follows. In Section II, we introduce a classification of the related work regarding time series forecasting. In Section III, we describe our proposed solution. Its performance evaluation and its comparison with existing approaches are the object of Section IV. The paper concludes in Section V.

II. Related Work

Several techniques have been proposed in the literature for function prediction, and more particularly for time series forecasting. They can be broadly classified into three main categories: statistical, machine learning-based and strategies mixing both types of methods. Although the list of techniques presented below is not exhaustive, this classification allows us to properly position the existing work and our approach.

A. Statistical-based approaches

Statistical approaches are traditionally used to make forecasts. One of the most widely recognized procedures for time series prediction is the Auto-Regressive Integrated Moving Average (ARIMA) one [START_REF] Yoo | Time-series forecast modeling on high-bandwidth network measurements[END_REF], which has been used in many studies. This technique has proven to be very effective when the data are stationary and show no upward or downward tendencies. However, the effectiveness of ARIMA remains relevant for linear systems only, making it inappropriate in many cases.

In [START_REF] Dai | EMD-Based Multi-Model Prediction for Network Traffic in Software-Defined Networks[END_REF], the Empirical Mode Decomposition (EMD)-based Multi-Model Prediction (EMD-MMP) algorithm has been proposed in order to predict network traffic in Software-Defined Networks. EMD-MMP combines the action of three methods: EMD [START_REF] Norden | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], [START_REF] Wu | Ensemble empirical mode decomposition: A noise-assisted data analysis method[END_REF], Auto-Regressive Moving Average (ARMA), and Support Vector Machine (SVM) -Regression (SVR) methods. Unlike ARIMA, this technique can be used to deal with non-linear signals and appears to be suitable for a variety of use cases.

Although statistical techniques are still relevant, they have recently been surpassed by machine learning-based ones [START_REF] Siami-Namini | A comparison of arima and lstm in forecasting time series[END_REF].

B. Machine learning-based approaches

Many recent studies are based on machine learning algorithms, which have the advantage of not requiring an explicit mathematical model. Moreover, once trained, the predictions are almost instantaneous, in contrast with statistical approaches. The LSTM neural network [START_REF] Hochreiter | Long Short-Term Memory[END_REF] is a top procedure in this category. It is based on a recurrent neural architecture that can process time series and is capable of learning and remembering input data over long sequences. This leads to its wide usage in time series prediction. In [START_REF] Siami-Namini | A comparison of arima and lstm in forecasting time series[END_REF], the authors compare the accuracy of ARIMA and LSTM for predicting time series data, resulting in a very significant superiority of LSTM over ARIMA.

The success of the LSTM algorithm in predicting or recognizing sequences has led to the emergence of several extensions. Among them, the Bidirectional LSTM algorithm is one of the most efficient variants [START_REF] Siami-Namini | A comparative analysis of forecasting financial time series using arima, lstm, and bilstm[END_REF]. The idea behind this technique is to consider not only the past values of the sequence but also the future ones, thus improving its prediction performance compared to a classical LSTM.

C. Combined approaches

Some approaches combine different techniques whether statistical or machine learning-based or even others. The M4 competition has shown that combined prediction methods have a higher accuracy than non-combined ones [START_REF] Makridakis | The M4 Competition: 100,000 time series and 61 forecasting methods[END_REF].

In [START_REF] Andreoletti | Network traffic prediction based on diffusion convolutional recurrent neural networks[END_REF], the authors proposed using a Diffusion Convolutional Recurrent Neural Network (DCRNN) to forecast network traffic. They used two layers of Diffusion Convolutional Gated Recurrent Unit (DCGRU). The comparison of this approach with an LSTM-based network, a CNN-based network, a CNN-LSTM-based network and a Fully Connected Neural Network show the superiority of DCRNN in forecasting congestion events.

The approach we propose in this paper falls into this category of methods, since we combine the FFT and an LSTM neural network.

III. A combined FFT/LSTM strategy for time-series prediction

In this section, we will introduce our approach, which consists in combining FFT data pre-processing with an LSTM neural network. In what follows, we detail the proposed architectures and the input/output conditioning to make prediction efficient, even for large measurement horizons.

A. Input/Output

Our work involves developing a component that receives old data as input to predict new data as output.

The data is a sequence of chronologically ordered values, each value d i representing the bandwidth read at time step i. From this data sequence, plus a sliding window d i+1 , d i+2 , . . . , d i+n of size n and the corresponding predicted value di+n+m (integer m ≥ 1 is called the prediction horizon), we generate a set of pairs (vector of values in sequence of size n, predicted value) where the vectors will be used as input to our component and the predicted value will be its output (see Fig. 1). We write

di+n+m = Φ(d i+1 , d i+2 , . . . , d i+n ) (1)
with Φ representing the prediction function.

We opted for a sufficiently large measurement horizon to avoid making multiple predictions before reaching the desired value, thus avoiding the propagation of errors. Indeed, given that we are at time step i + n, predicting the value at time i+n+m implies the prediction of all the values before this last one. This implies a prediction from inputs with successive estimation errors. To overcome this problem, we proposed to estimate the value at time step i + n + m directly, without resorting to the intermediate values.

d 1 , d 2 , . . . , d i , d i+1 , d i+2 , . . . , d i+n , d i+n+1 , . . . , d i+n+m , . . .

Sliding window

Predicted value Fig. 1: Sliding window and predicted value.

B. Architectures

In our approach, we have combined LSTM and FFT (Fast Fourier Transform) in different ways. In what follows, we will present the DFT and the FFT and the different ways they have been used in our architectures.

1) DFT and FFT: The Discrete Fourier Transform (DFT) is a function allowing to perform a Fourier analysis on digitized signals. The DFT of a signal [x 0 , x 1 , ..., x N -1 ] of length N , is another signal [y 0 , y 1 , ..., y N -1 ] defined by

y k = N -1 ∑ n=0 x n [ cos ( 2π N kn ) -i sin ( 2π N kn )] (2)
where k ∈ {0, 1, ..., N -1} and i is the imaginary unit. The DFT of a vector is calculated by means of the Fast Fourier Transform (FFT) [START_REF] Cooley | An Algorithm for the Machine Calculation of Complex Fourier Series[END_REF], which is a fast algorithm with complexity O(n log n) for a vector of size n. The DFT is used in several domains, such as digital treatment of signals, image compression, voice recognition, etc. As for the classic Fourier transform, it converts signals from the time domain into the frequency domain. Note that the transformation to the frequency domain with DFT causes the loss of temporal information, but prediction needs this temporal dimension. We can solve this problem by using a sliding window. Indeed, the DFT provides a frequency dimension and the sliding window allows keeping the time dimension.

2) Architectures: We have produced four different architectures in our studies: LSTM (method 1), RealFFT-LSTM (method 2), CombinedFFT-LSTM (method 3) and HSFFT-LSTM (method 4), with two versions each, a simple LSTM layer and a double-layer LSTM, see Fig. 2. In the first one, the sliding windows are passed directly to an LSTM neural network to predict the future value. The sliding window, in the second architecture, undergoes a DFT preprocessing, which returns a vector of complex numbers, and then their real parts are injected into the LSTM input. In the third architecture, real and imaginary parts of the DFT vector are concatenated before being injected to the LSTM input. In the fourth proposed architecture, we have combined the HSFFT [13] 1 , which is explained below, and an LSTM.

The HSFFT module consists in duplicating the signal in order to make it symmetrical. Thus, we have to consider that the sliding window represents only half of the input signal and we duplicate it in order to obtain a symmetrical signal as shown in Fig. 3. This balanced signal is processed by the FFT algorithm, which provides a true symmetrical frequency output vector having a zero imaginary part. Then, only a part of this frequency vector is transmitted to the LSTM (duplicated information is removed).

The choice of HSFFT over FFT was made after conducting some comparative tests between the two methods. The results showed that HSFFT reduces the error rate, especially when combining two LSTMs.

IV. Performance evaluation A. Simulations' settings

In the following, we compare the different versions of the LSTM method that we have previously introduced. We considered the following methods: the LSTM, the RealFFT-LSTM, the CombinedFFT-LSTM and the HSFFT-LSTM. They are composed of one or two LSTM layers with the same number of neurons. A fifth module has been devoted to the ARIMA method. All modules were developed in Python using the Keras library [5] with a Tensorflow backend [START_REF] Tensorflow | TensorFlow: Large-scale machine learning on heterogeneous systems[END_REF]. Frequencies that have a zero imaginary part (vector of size 2n-1)

Fig. 3: Detail of the HSFFT module.

We performed the tests on real network traffic data coming from a mobile network operator. The benchmark contains 8928 values of bandwidth measurement at 10-minute intervals, see Fig. 4. Data was first normalized to be compliant with the LSTM input. We have considered a sliding window with a size varying from 20 to 70 values (Fig. 1).

Data was divided into two parts: 70% was used for the learning phase and 30% for the test phase. For the training, we set the batch size to 1 and to learn, we use an Adam optimizer with a learning rate of 0.01. To evaluate the effectiveness of each module, we considered the calculation of the MAPE (Mean Absolute Percentage Error) between the real values and the predicted ones. The comparative study was conducted in two phases:

1) Phase 1: It consisted of comparing the HSFFT-LSTM method with similar ones, using a deep learning approach (LSTM, CombinedFFT-LSTM and RealFFT-LSTM). 2) Phase 2: Here, we compared our approach to the ARIMA method.

The obtained results of each module are treated, analyzed, and represented by box plots. We performed 72 series of tests in which every test was executed 20 times (for the calculation of confidence intervals). The tests have been carried out by changing the most influencing parameters one at a time in order to determine the error rate for different prediction horizons, and using different values for the size of the sliding window.

We have taken into consideration the effectiveness of each module based on the following criteria:

• accuracy over the prediction horizon (by varying the number of neurons within the LSTM and the size of the sliding window);

• accuracy as a function of the size of the sliding window;

• accuracy as a function of the number of LSTM neurons;

• prediction latency.

B. Results 1) Accuracy over the prediction horizon: We performed a first test comparing the LSTM, RealFFT-LSTM, CombinedFFT-LSTM and HSFFT-LSTM methods according to the prediction horizon, for an LSTM of 32 neurons, a sliding window of size 30 and a prediction horizon of 1 hour, 2 hours and 3 hours. We used the MAPE error, plotted in Fig. 5. The obtained results show that increasing the prediction horizon has an almost linear impact on the prediction error. The increase of this error with the prediction horizon depends, however, on the strategy applied. Indeed, the larger the prediction horizon, the larger the difference between the considered strategies. Besides, one can notice that the methods using filters (i.e., RealFFT-LSTM, CombinedFFT-LSTM and HSFFT-LSTM) clearly outperform the LSTM alone, when increasing the distance with the forecasted value. For a prediction with a distance of 1h (6th value in the future), the different approaches show, nevertheless, results with the same order of magnitude. The two methods CombinedFFT-LSTM and HSFFT-LSTM exhibit the best performance with a slight advantage to the second one.

2) Accuracy as a function of the size of the sliding window: The first tests led us to conduct another series of experiments to analyze the behavior of the LSTM and HSFFT-LSTM methods when we changed the size of the sliding window. We conducted 6 series of 20 tests each, for different values of the sliding window, whose results are shown in Fig. 6. We also performed these experiments with two combined LSTM layers (see Fig. 7).

When using only one LSTM layer, we notice that for a rather small window size, HSFFT-LSTM shows a better accuracy compared to the LSTM alone. However, when the sliding window exceeds a certain size both methods give roughly the same performance with a slight advantage to 3) Accuracy as a function of the LSTM neurons' number: We analyze, here, a series of tests designed to examine the impact of the number of neurons within the LSTM network on the error rate, over multiple measurement horizons. The comparison is made only between the two methods LSTM and HSFFT-LSTM. We chose an LSTM with 32, 64, 96 and finally 128 neurons, horizons of one hour, two hours and three hours, and each time we made 20 tests (see Fig. 8).

The results of this series led us to the conclusion that HSFFT-LSTM is better than an LSTM alone, regardless of the number of neurons. Moreover, with only a few neurons, HSFFT-LSTM still makes excellent predictions, unlike the LSTM, which deteriorates significantly for horizons far in the future. The HSFFT-LSTM procedure remains powerful for a fairly large prediction horizon, and even with a reduced number of neurons, which allows saving system resources.

4) HSFFT-LSTM compared to ARIMA with different sliding windows: We compare, here, our approach with a classical statistical method, which led us to perform another series of prediction tests using the ARIMA and the HSFFT-LSTM modules. Our sliding window size was 20, 30 and 40, using prediction horizons of one hour, two hours and three hours.

The results of these tests conducted us to deduce that our approach is especially good for large prediction horizons and large sliding windows (see Fig. 9).

5) Prediction latency:

To compare the execution time of the ARIMA, LSTM and HSFFT-LSTM techniques, we conducted some experiments with a sliding window of size 10, and we executed ARIMA using data of size 288 to ensure its convergence. We divided the ARIMA execution time by 30 to make a fair comparison with 10-sized sliding window HSFFT-LSTM.

We found that LSTM and HSFFT-LSTM are significantly better in terms of execution time compared to In conclusion, we have seen that the inclusion of the FFT allows a significant improvement in the prediction of values over various measurement horizons. Even if it is difficult to rigorously justify the underlying reasons, we believe that the shift from the time domain to the frequency domain, with amplitudes reflecting the importance of the frequencies, allows the LSTM to better grasp the parameters impacting the prediction the most and to separate them more easily from potential noises that are not very visible on the time scale.

V. Conclusion

In this work we propose the use of data preprocessing before injecting them into an LSTM neural network for time series prediction. Our approach is based on the use of an FFT filter applied to a sliding window of data.

We combined LSTM with FFT in several ways and performed a series of comparative tests, exploring the procedures' performances by changing several parameters. Our study led us to conclude that FFT-LSTM is more effective in terms of prediction horizon and execution time than the classical ARIMA. Compared to similar methods such as LSTM alone, our analysis showed that for a large dataset and with a reduced number of neurons, our approach gives good results even for a large prediction horizon. In other words, the combination that we propose of methods not only reduces significantly the prediction error, but can also reduce the number of required neurons, thus accelerating the learning process. In addition, our approach maintains good performance over long prediction horizons.

In the future, we intend to use this technique in the management of an SDN/NFV network, notably through the dynamic scaling of services. We also plan to improve its rapid adaptation to potential variations in traffic patterns.

Fig. 4 :

 4 Fig. 4: Network data traffic.

Fig. 5 :

 5 Fig. 5: Comparison between LSTM-based methods according to different prediction horizons (sliding window size = 30, neurons number = 32).

Fig. 6 :Fig. 7 :

 67 Fig. 6: Comparison according to different sliding window size (horizon = 3 hours) when using 1 LSTM layer.

Fig. 8 :

 8 Fig. 8: Comparison between LSTM and HSFFT-LSTM methods according to the number of neurons and different prediction horizons (sliding window = 30).

Fig. 9 :Fig. 10 :

 910 Fig. 9: Comparison between ARIMA and HSFFT-LSTM methods. .

University Oran 1 RIIR Laboratory, Oran, Algeria

Univ. Rennes, Inria, CNRS, IRISA

Inria, Univ. Rennes, CNRS, IRISA


Diego Abalos
email: d.abalos@agro.au.dk

Tatiana F Rittl

Sylvie Recous

Pascal Thiébeau

Cairistiona Frances

Kees Elizabeth Topp

Jan Van Groenigen

Klaus Butterbach-Bahl

Rachel E Thorman

Kate E Smith

Ishita Ahuja

Cairistiona F E Topp

Jan Kees

Van Groenigen

Jørgen E Olesen

Marina A Bleken

Robert M Rees

Sissel Hansen

Elena Paoletti

Predicting field N2O emissions from crop residues based on their biochemical composition: A meta-analytical approach

Keywords:

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Crop residue production exceeds five billion tons per year worldwide [START_REF] Cherubin | Crop residue harvest for bioenergy production and its implications on soil functioning and plant growth: a review[END_REF]. A large part of this material is returned to agricultural land as a strategy to restore soil organic matter stocks and build soil fertility. Indeed, together with the return of animal excreta, the incorporation of crop residues is the most important method for farmers to maintain and sustain soil organic matter [START_REF] Powlson | Carbon sequestration in european soils through straw incorporation: limitations and alternatives[END_REF]. Although returning crop residues has many benefits for soil health, it may also increase the emission of the powerful greenhouse gas nitrous oxide (N 2 O) (IPCC, 2006;[START_REF] Chen | Increased N2O emissions during soil drying after waterlogging and spring thaw in a record wet year[END_REF][START_REF] Hansen | Reviews and syntheses: review of causes and sources of N2O emissions and NO3 leaching from organic arable crop rotations[END_REF]. Nitrous oxide emissions from agricultural land account for approximately 60% of the total global anthropogenic N 2 O emissions [START_REF] Reay | Global agriculture and nitrous oxide emissions[END_REF], and are therefore considered a key target in the Paris Agreement in order to limit global warming to 1.5 °C or 2 °C [START_REF] Clark | Global food system emissions could preclude achieving the 1.5°and 2°C climate change targets[END_REF].

Although crop residues generally increase N 2 O emissions, negative and no effects have also been reported [START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF]. The effect of residue retention on N 2 O is difficult to predict since it is mediated via multiple processes. Crop residues might increase N 2 O emissions by providing N for soil microbial processes [START_REF] Li | Carbon sequestration in arable soils is likely to increase nitrous oxide emissions, offsetting reductions in climate radiative forcing[END_REF], including nitrification and also denitrification [START_REF] Surey | Potential denitrification stimulated by water-soluble organic carbon from plant residues during initial decomposition[END_REF]. Accordingly, N 2 O emissions from crop residues are often predicted using a fixed emission factor based on the amount of nitrogen (N) returned to the soil with the residues (IPCC, 2019). However, this approach neglects the impact of residue-C addition on microbial dynamics and denitrification processes. Residue C:N ratio, a commonly used index for residue quality, accounts for net N immobilization or mineralization from crop residues, and is also used as predictor of residue-derived N 2 O emissions [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF]. Negative relationships between residue C:N ratio and N 2 O emissions are frequently reported [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF], and are generally attributed to decreased N availability due to microbial immobilization of soil mineral N with high C:N residues [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF]. However, this could be an oversimplification, because crop residues differing in their C:N ratio are likely to differ also in terms of soluble and structural contents, which determine the extent and dynamics of residue decomposition and consequently also affect N 2 O production. For example, increases in lignin and polyphenols concentration can be negatively related to N 2 O emissions [START_REF] Millar | Chemical composition, or quality, of agroforestry residues influences N2O emissions after their addition to soil[END_REF][START_REF] Garcia-Ruiz | N2O emission from soil following combined application of fertiliser-N and ground weed residues[END_REF], whereas water-extractable organic C may stimulate denitrification [START_REF] Surey | Potential denitrification stimulated by water-soluble organic carbon from plant residues during initial decomposition[END_REF]. The direct links between biochemical crop residue characteristics and N 2 O emissions, and particularly the correlations between them, remain poorly understood.

Predicting N 2 O emissions from crop residues is also challenging because these emissions depend on multiple interactions between crop type, environmental conditions, and soil properties [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Lehtinen | Effect of crop residue incorporation on soil organic carbon and greenhouse gas emissions in european agricultural soils[END_REF][START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF]. For example, N 2 O emissions from crop residues are negatively and linearly related to soil clay content [START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF], because high clay content is associated with lower soil porosity and aeration, and thus lower oxygen availability [START_REF] Skiba | Emissions of nitric oxide and nitrous oxide[END_REF]. However, N 2 O emissions from crop residues are not linearly related to soil pH and SOC, instead peaking at intermediate values of either of those soil properties [START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF]. Residue return stimulates N 2 O emissions more strongly in regions with a warmer climate (e.g., warm temperate and subtropical) compared to cool temperate zones, probably due to increased decomposition rate. Residue-derived N 2 O emissions may be negatively associated with mean annual precipitation [START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF]; yet, other studies showed increased emissions at intermediate soil water content (i.e. water filled pore space 60-90%; [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF], perhaps due to confounding effects via soil type. While previous meta-analyses have evaluated these factors individually, their hierarchical importance has been overlooked.

The objective of this meta-analysis was to synthesise the results of field studies measuring the effect of crop residue retention on soil N 2 O emissions. We aimed to identify the main drivers of soil N 2 O emissions associated with the application of crop residues and to rank their importance. Our study differs from previous meta-analyses in several aspects: we only collected information from field studies to avoid unrealistic emissions from crop residues; used the most comprehensive approach so far to gather biochemical characteristics of crop residues; used a random-meta-forest approach, which accounts for multiple drivers simultaneously including nonlinear relationships; and developed and tested for the first time a simplified classification based on the stage of residue maturity. We hypothesized that the relationships between crop residue biochemical characteristics and N 2 O emissions can be unfolded to develop a predictive framework of crop residue effects on N 2 O emissions.

Materials and methods

Data collection and data gap filling

Crop residue and soil N 2 O emissions database

We collected data on N 2 O emissions from field experiments with crop residue retention treatments across the world. We defined crop residues as any constituent of a crop (above-and below-ground) returned to the soil. These may originate from cash crops, cover crops, green manures, grassland or others. Depending on the crop, the residues may include below-ground and above-ground biomass, except for root crops for which it refers to above-ground biomass only since most of the roots are harvested. The crop residues treated through biogas fermentation, composting and other methods were not included. Another difference with respect to previous meta-analyses is that we included grassland termination data in our database, because temporary grasslands make up a substantial share of arable land in Europe, especially in the Scandinavian countries and in some regions with high livestock densities, e.g. the Netherlands, Belgium and Bretagne, France [START_REF] Lesschen | Task 1defining and classifying grasslands in Europe[END_REF]. Peer-reviewed research articles published before 1 January 2021 were searched in Google Scholar and Web of Science with the following keywords: "residue" OR "straw" OR "stubble" AND "nitrous oxide" OR "N 2 O" OR "greenhouse gas" OR "GHG" OR "emission factor". The initial search resulted in 337 publications. We also included results from two unpublished studies conducted by authors of this meta-analysis (Bleken et al., n.d.; Ernfors et al., n.d.). All publications were screened and included in the database if they met the following criteria: (i) the experimental design had to be sufficiently detailed regarding the treatments and recent history; (ii) included treatment replicates; and (iii) were conducted under realistic field conditions. All the studies included pairwise comparisons in which all factors (e.g., fertilization, irrigation, tillage practices) were the same for the treatment (with crop residue retention) and control (with crop residue removal) groups. For grassland, we compared permanent grassland (control) to grassland termination (treatment), the analogous practice to crop residue incorporation in these systems. Note that the control treatments do not always mean residue-free soil and its definition varies among studies, e.g. in green manure studies, the control may be a treatment without the addition of herbage with the presence of roots or a fallow plot where the roots were removed. The selected studies (n = 78) provided 367 pairwise comparisons. Cumulative N 2 O emissions (g N 2 O-N ha -1 ) for each crop residue treatment and control were extracted together with a measure of variance and the number of replicates. The length of the studies ranged from 5 to 575 days. To account for potential biases due to the inclusion of shortterm studies (<40 days), we ran an initial sensitivity analysis excluding short-term studies and it did not affect our main results. For each study in our dataset, we tabulated details on the site and experimental conditions. When data were only presented graphically, WebPlotDigitizer was used to extract data points. If relevant information was not reported, the authors were contacted to supply the missing information.

Data aggregation

We identified 21 factors as potential drivers of N 2 O emissions from crop residues (Table 1). These factors related to crop type, residue type, residue properties, soil properties, and climatic conditions, and included both categorical and continuous factors. A detailed data overview is presented in a publicly available database accompanying this publication (Rittl et al., n. d.). For crop type, we distinguished between cereal (mature harvested cereals and rapeseed), rice, cover crop (including immature cereals, grasses, clovers, Brassica spp. cultivated as cover crops or one-year green manures), grassland, grain legumes, vegetables, and sugar cane, and double cropping (combination of two species at different biomass ratios). For soil properties, soil texture was divided into clay (clay, sandy clay, silt clay), loam (loam, clay loam, silty clay loam, silt loam, sandy loam, sandy clay loam, fine silt, silt), and sand (sand, loamy sand, fine sand) [START_REF] Quemada | Metaanalysis of strategies to control nitrate leaching in irrigated agricultural systems and their effects on crop yield[END_REF]. Soil pH values measured in CaCl 2 were converted to H 2 O-based values using the equation proposed by [START_REF] Ros | Towards optimal use of phosphorus fertiliser[END_REF]. Soil pH values in our dataset ranged from 4.3 to 9.3, and were grouped into acidic (<6), neutral (6-7), and alkaline (>7) soils. We converted organic matter content to organic C content using the standard conversion factor of 0.58 [START_REF] Pribyl | A critical review of the conventional SOC to SOM conversion factor[END_REF].

We calculated the Aridity Index (AI) as the ratio of mean annual precipitation to mean annual potential evapotranspiration (UNEP). The long-term mean annual precipitation was extracted from the publications, and the long-term mean annual potential evapotranspiration for each site was extracted from the Global Aridity Index geodatabase (https://cgiarcsi. community/data/global-aridity-and-pet-database/). We used the classification of wet (AI >1) and dry land (AI <1) utilized by IPCC (2019) to differentiate N 2 O emission factors from N sources (including organic amendments such as crop residues).

We assessed the impact of study duration on a thermal time basis, using the normalized time approach proposed by [START_REF] Andren | Barley straw decomposition in the field: a comparison of models[END_REF], and applied successfully to soil C and N mineralization, gross N fluxes and C-N models [START_REF] Aita | Short-term kinetics of residual wheat straw C and N under field conditions: characterization by 13C15N tracing and soil particle size fractionation[END_REF][START_REF] Brisson | Conceptual Basis, Formalizations and Parameterization of the STICS Crop Model Versailles[END_REF]. In short, the total duration of the experiment was converted to number of days at a reference temperature of 15 °C (i.e. ND15), following the approach described by [START_REF] Thiébeau | Dynamiques de décomposition des résidus de cultures sur des exploitations pratiquant l'agriculture de conservation en région Grand Est, France[END_REF].

Crop residue quality characteristics and categorical groups

From each study in our dataset, we extracted information on the chemical characteristics of crop residues. We calculated the N returned with crop residues from the residue N concentration and the dry matter (DM) added. When C and N added were reported, we calculated crop residue C:N ratio. When these data were not provided, mean values were taken from the crop residue quality dataset provided by [START_REF] Thiébeau | Dataset of biomass and chemical quality of crop residues from european areas[END_REF] for the same species or group of species. Residue concentration of soluble fraction (neutral detergent soluble, referred to as soluble NDS fraction), hemicellulose, cellulose, and lignin (according to Van Soest (1963)) were tabulated when these were reported. When these data were not available, which was the case for the vast majority of studies included in the database (>90%), we used crop residue quality data from [START_REF] Thiébeau | Dataset of biomass and chemical quality of crop residues from european areas[END_REF]. Previous research has shown that, on a global scale, most of the variation in plant characteristics is represented by species identity [START_REF] Kattge | TRYa global database of plant traits[END_REF]. Accordingly, the use of a global dataset for residue biochemical characteristics was deemed appropriate for our global meta-analysis. Hemicellulose, cellulose and lignin relative contents (% total DM), indicative of the composition of the insoluble residue fraction (i.e. the plant cell walls), were used, the difference between their sum and 100%, being the soluble NDS fraction. The lignocellulose index (LCI) [lignin: (lignin + cellulose + hemicellulose)] was used as a criteria to express the recalcitrance of the plant cell wall [START_REF] Herman | The relationship between rates of lignin and cellulose decay in aboveground forest litter[END_REF]. Water soluble C (WSC) determined after water extraction, was expressed as a percentage of total C. Crop C:N ratio was divided into three categories (<20, 20-60, >60).

We used two criteria to categorize crop residues: according to maturity stage, and to the type of residue generated. Crop residues were classified as mature or immature, based on the stage of physiological maturity of the plant at which the residue were generated [START_REF] Sylvester-Bradley | Minimising nitrous oxide intensities of arable crop products (MIN-NO)[END_REF]. This could be by cultivation practice (harvesting of annual crops at the end of the cycle, harvesting of root or vegetable crops, mechanical destruction, grassland mowing or grassland renewal) or possibly naturally (senescence), irrespective of the plant type. Incorporation of immature residues into the soil happens during: (i) destruction of green cover crops or [START_REF] Sylvester-Bradley | Minimising nitrous oxide intensities of arable crop products (MIN-NO)[END_REF]. In our database, the classification of immature residues refers to cover crops, vegetable residues and grassland renewal. Residues were classified as mature for a range of post-crop harvest residues such as straw and stem residues from cereals, rice, grain legumes, oilseed, sugar cane and others. Second, we categorized studies according to the type of residue generated (Table 1). Residues types were grouped into mature above-ground biomass, straw, green plant biomass, and senescent plant biomass (Table 1). Mature above-ground biomass encompassed above-ground residues at maturity for all crops except cereals. Straw indicates above-ground residues of cereals including rice at maturity. Green plant biomass is the total plant harvested at vegetative stage (i.e. cover crops, grassland renewal, vegetables). Senescent plant biomass refers to leaves of perennial crops, in our case only sugar cane.

Dataset syntheses

The database included studies from 69 sites in 20 countries. Europe contributed 41% of the comparisons (n = 148). For Europe, 32% of the data were from residues classified as cover crop; 80% as green plant biomass; 60% of the N 2 O comparisons were from experiments where fertilizer was applied to either the growing crop generating the residue, or at some point following residue return. Regarding climate, 63% of our comparisons originated from field sites in dry climates (AI <1), and 37% from sites in wet climates (AI >1). The range of each continuous factor is shown in Table 1.

Data analysis

We used the log response ratio (LnRR) as effect size, which is a common metric in meta-analyses [START_REF] Hedges | The meta-analysis of response ratios in experimental ecology[END_REF][START_REF] Osenberg | Resolving ecological questions through meta-analysis: goals, metrics, and models[END_REF]. We performed a weighted mixed-effects meta-analysis, using the rma.mv function of the 'metafor' package [START_REF] Viechtbauer | Conducting meta-analyses in R with the metafor[END_REF], including "study" as a random effect because several studies contributed more than one effect size. Effect sizes from individual comparisons were weighted by the inverse of their variance. Missing variances were estimated using the average coefficient of variation across the dataset [START_REF] Van Groenigen | Sequestering soil organic carbon: a nitrogen dilemma[END_REF]. We used a Wald test to evaluate statistical differences between subgroups within categories at p-value < 0.05.

We used a random-meta-forest approach to identify the most important predictors of crop residue effects on N 2 O emissions. This approach is based on the machine learning random forest algorithm. This approach incorporates the variance and weight of each study as in classic meta-analysis, is robust to overfitting, allows for numerous predictors, and accounts for nonlinear relationships [START_REF] Terrer | Nitrogen and phosphorus constrain the CO2 fertilization of global plant biomass[END_REF][START_REF] Terrer | A trade-off between plant and soil carbon storage under elevated CO2[END_REF]). Among the 21 potential factors, we discarded those based on expert knowledge to avoid possible artifacts associated with arbitrary category definitions (i.e. crop type, residue type, maturity class) (Table 1), and those with more than 30% missing data (soil texture, clay, annual temperature, soil bulk density, soil total N, and ND15). The final dataset consisted of 267 observations with information on 12 predictors, which are indicated in Table 1. Using 'metaforest' (Van Lissa, 2017), we performed variable pre-selection with 10,000 iterations and 100 replications, using a recursive algorithm in the preselect function from metafor. Predictors with negative variable importance were dropped using the preselect_vars function. Predictors with positive predictive performance were then used to optimize the model. Parameters of the meta-forest model were optimized using the train function from the caret package [START_REF] Kuhn | Building predictive models in R using the caret package[END_REF], and calculated 10-fold cross validated R 2 with 75% of the data used as training data and 25% as validation data. We generated partial dependence plots to depict the association of each selected predictor with the effect size, while accounting for the average effect of all other predictors.

To check the robustness of our random-meta-forest approach, we ran another mixed-effects meta-regression model selection procedure with the six most important predictors identified by random forest analysis. The model selection was based on maximum likelihood estimation, which was carried out with the 'glmulti' package in R [START_REF] Calcagno | Glmulti: an R package for easy automated model selection with (generalized) linear models[END_REF]. The importance of each predictor was computed as the sum of Akaike weights for models that included this predictor. A cutoff of 0.8 was set to differentiate between essential and non-essential predictors [START_REF] Chen | Long-term nitrogen loading alleviates phosphorus limitation in terrestrial ecosystems[END_REF][START_REF] Terrer | Nitrogen and phosphorus constrain the CO2 fertilization of global plant biomass[END_REF].

Nitrogen fertilizer application is one of the most important drivers of N 2 O emissions, and it may affect the relationships between our identified controlling factors and N 2 O emissions from crop residues, as observed for organic amendments [START_REF] Charles | Global nitrous oxide emission factors from agricultural soils after addition of organic amendments: a meta-analysis[END_REF] and for crop residues [START_REF] Shan | Effects of crop residue returning on nitrous oxide emissions in agricultural soils[END_REF]. Our initial analyses confirmed the importance of N fertilizer addition (p < 0.01). Thus, we used 'metafor' to test for interactions between fertilizer application (i.e. with or without application) and each of the most important predictors identified by the random-meta-forest approach. When the interaction term was non-significant, we proceeded with the analyses using the complete dataset; when the interaction was significant, we conducted separate analyses for studies with and without fertilizer application. We fitted linear and quadratic meta-regressions to identify the best model describing the relationships between the most important predictors and LnRR, based on the Akaike information criterion (AIC); the model with lower AIC was retained.

Principal component analyses (PCA) were performed with R (prcomp available in the base package, and the 'ggbiplot' package; [START_REF] Vu | ggbiplot: a ggplot2 based biplot[END_REF] to visualize the main biochemical characteristics of mature and immature residues, and of the main crop type.

Results

Overall effect of crop residues on N 2 O emissions

Averaged across all studies, crop residue retention increased N 2 O emissions by 43% (95% CI: 23-66%; p < 0.001). However, the differences between crop types were large (Fig. 1). Crop residues increased N 2 O emissions for cover crops (50%, CI: 19-88%), grasslands (216%, CI: 114-366%) and vegetables (183%, CI: 85-334%), but not for cereals, grain legumes, rice or sugar cane (Fig. 1).

Main predictors of crop residue induced N 2 O emissions

Our random-meta-forest approach identified N returned with crop residues, WSC, cellulose, lignin, hemicellulose and soluble NDS fraction as the most important predictors of crop residue effects on N 2 O emissions (Supp. Fig. 1). Model selection analysis based on Akaike weights confirmed the importance of these predictors, as all of them reached the threshold of 0.8 (Supp. Fig. 2). The importance of residue quality properties was therefore higher than that of the soil properties and climatic factors included in the random-meta-forest analysis (Table 1).

The amount of N returned (kg N ha -1 yr -1 ) with the crop residues affected the N 2 O emissions (AIC = 481.49, n = 298), and the overall effect of the amount of residue-N was not modified by the application of fertilizer N (p > 0.05). Nitrous oxide emissions from crop residues increased with the amount of residue-N returned, reaching a maximum at ~200-250 kg N ha -1 yr -1 applied (Fig. 2). Crop residue-derived N 2 O emissions decreased with increasing values of cellulose, whereas intermediate hemicellulose values (~30%) were associated with the lowest LnRR (Fig. 2). The effects of cellulose (AIC = 592.87, n = 346) and hemicellulose (AIC = 612.05, n = 345) on N 2 O emissions from crop residues were not affected by fertilizer application either.

Nitrogen fertilizer application affected the relationship between crop residue N 2 O emissions and WSC (with fertilizer: AIC = 208.02, n = 181; without fertilizer: AIC = 133.58, n = 70), soluble NDS fraction (with: AIC = 297.59, n = 210; without: AIC = 148.13, n = 79), and lignin (with: AIC = 313.18, n = 213; without: AIC = 153.08, n = 84), but the trends were in general consistent with and without fertilization (Fig. 3). N 2 O emissions from crop residues increased with WSC values and when the residue soluble NDS fraction exceeded 40%. High lignin content in the crop residues corresponded with smaller treatment effects on N 2 O emissions for studies without added fertilizer, but not for studies with added fertilizer.

Residue type and N 2 O emission

Maturity class had a strong influence on the magnitude of N 2 O emissions (p < 0.001), with greater emissions from immature than from mature residues (Wald-type test, p < 0.001) (Fig. 4). The type of residue significantly affected crop residue derived N 2 O emissions (p < 0.001). Green plant biomass increased N 2 O emissions by 89%, and mature aboveground biomass by 40% (Supp. Fig. 3). Conversely, straw and senescent plant biomass did not increase N 2 O emissions.

Predictors with lower importance

The response of N 2 O emissions to crop residue incorporation was not explained by the lignocellulose index (LCI), annual mean precipitation (i.e. these two predictors were not selected by the random-meta-forest approach), soil texture, soil bulk density, annual mean temperature or study length either expressed by calendar day or by standard day (ND15) (i.e. these predictors had no significant effect on LnRR based on mixed-effects meta-analysis). However, simple linear meta-regressions indicated that residue C:N ratio (p < 0.001), soil pH (p = 0.014), SOC (p = 0.011), soil N (p = 0.002) and Aridity Index (p = 0.021) exerted a small but significant influence on crop residue N 2 O emissions (Supp. Figs. 4,5). High C:N ratios (above c. 60), neutral soil pH, low SOC and low soil N content led to lower N 2 O emissions. The effect size of applying crop residues in humid climates was 2 times larger than in dry climates.

Relationship between residue biochemical characteristics

A PCA analysis was used to determine the association between residue biochemical characteristics, grouped for maturity class and crop type (Fig. 5). The first two principal components accounted for 72.4% of the total variance. We found strong relationships between residue biochemical characteristics: WSC, NDS fraction and residue N returned were positively correlated among them, and negatively related to cellulose, hemicellulose and C:N ratio. Immature residues were characterized by high WSC, NDS fraction and N returned with low cellulose, hemicellulose, lignin fractions, LCI and C:N ratio; the opposite was true for mature residues (Fig. 5a). Cover crop residues generally contained a high soluble NDS fraction and WSC, vegetable residues returned relatively large amounts of N, and grassland residues presented low lignin concentration (Fig. 5b). Cereals and sugar cane residues were characterized by a high cellulose and hemicellulose content, and a high C:N ratio.

Discussion

Magnitude of N 2 O emissions as affected by crop type

In line with other meta-analyses [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF][START_REF] Muhammad | Regulation of soil CO2 and N2O emissions by cover crops: a meta-analysis[END_REF][START_REF] Shan | Effects of crop residue returning on nitrous oxide emissions in agricultural soils[END_REF], we found a significant stimulation of soil N 2 O emission following the application of residues from vegetables, cover-crops and grasslands renewal, but not from cereals, grain legumes, rice and sugar cane. For vegetables, we observed an average increase in N 2 O emission of around 183%, which is higher than reported in the meta-analysis of [START_REF] Shan | Effects of crop residue returning on nitrous oxide emissions in agricultural soils[END_REF] for lettuce and bean residues (123% and 138% increase, respectively). Our finding that incorporation of cover crop residue stimulated N 2 O emissions is consistent with [START_REF] Muhammad | Regulation of soil CO2 and N2O emissions by cover crops: a meta-analysis[END_REF], although they only found increased N 2 O emissions for legume-based cover crops. As we did not distinguish between residues from cover crops with and without legumes, we cannot directly compare our results with theirs.

The main drivers behind crop residue N 2 O emissions

Residue biochemical characteristics were strong predictors of N 2 O emissions from crop residues. Easily degradable fractions (WSC, soluble NDS fraction), structural fractions (cellulose, hemicellulose and lignin), and N returned with crop residues (resulting from residue N concentration and biomass) had the most predictive power. Crop residues rich in easily degradable fractions increased N 2 O emissions, probably by providing C sources for biological N transformation processes and by accelerating microbial growth [START_REF] Sahrawat | Nitrous oxide emission from soils[END_REF]. The easily degradable fractions can also promote N 2 O emissions by creating an anaerobic environment, and by providing energy sources for nitrifiers, and electron acceptors for bacterial and fungal denitrifiers [START_REF] Palatinszky | Cyanate as an energy source for nitrifiers[END_REF][START_REF] Surey | Potential denitrification stimulated by water-soluble organic carbon from plant residues during initial decomposition[END_REF]. Indeed, during short-term anoxic conditions due to high O 2 demand from microbial activity, denitrifying organisms rely on soluble compounds from fresh residues as major C sources [START_REF] Surey | Potential denitrification stimulated by water-soluble organic carbon from plant residues during initial decomposition[END_REF]. Conversely, structural fractions, which are relatively resistant to decomposition, were negatively related to N 2 O emissions from crop residues. This is consistent with a large body of literature showing that recalcitrant materials such as stubbles and roots can take longer to mineralize, and can lead to net immobilization of soil mineral N (e.g. [START_REF] Kuzyakov | Priming-effekte durch rhizodeposition von pflanzen[END_REF][START_REF] Gentile | Interactive effects from combining fertilizer and organic residue inputs on nitrogen transformations[END_REF]. The stimulatory effect of N returned with crop residues on N 2 O emissions is well established [START_REF] Li | Evidence for denitrification as main source of N2O emission from residue-amended soil[END_REF][START_REF] Hansen | Reviews and syntheses: review of causes and sources of N2O emissions and NO3 leaching from organic arable crop rotations[END_REF], as the release of residue N after mineralization increases the availability of soil mineral N for nitrification and denitrification, which are likely to be source processes under different environmental conditions.

The divergent relationships between individual biochemical characteristics and N 2 O emissions makes predicting N 2 O emissions from crop residues challenging. However, our results shed light on how to overcome these difficulties. This is because crop residue biochemical characteristics are correlated, as shown by the multivariate analysis (Fig. 5). From a quantitative point of view, this is of course linked to the characterization used (i.e. Van Soest method), with values representing relative proportions in fractions of increasing recalcitrance, expressed as a function of the total mass [START_REF] Soest | Use of detergents in the analysis of fibrous feeds. II. A rapid method for the determination of fiber and lignin[END_REF][START_REF] Soest | Use of detergents in the analysis of fibrous feeds. IV. Determination of plant Cell-Wall constituents[END_REF]. The biological explanation is that the distribution of the different pools within the plant tissues changes with plant development, with characteristic compositions linked to given developmental stages and functions (e.g. leaf vs. stem vs. root). For example, residues from cereals harvested at mature stages have less soluble and more structural fractions than cereals harvested earlier in the growth cycle [START_REF] Bertrand | Soil decomposition of wheat internodes of different maturity stages: relative impact of the soluble and structural fractions[END_REF], causing a negative association between soluble and structural fractions (Supp. Fig. 6). Our results also show a strong positive correlation between residue N returned and the soluble NDS fraction, confirming previous work [START_REF] Schmatz | How the mass and quality of wheat and vetch mulches affect drivers of soil N2O emissions[END_REF]. These biochemical associations indicate that crop residues can be categorized beyond single residue properties (e.g. N content, cellulose, hemicellulose, lignin), and that aggregate indicators of residue composition are probably the most useful predictors of residue effects on soil N 2 O emissions.

Maturity class: a new and robust predictor of N 2 O emissions from crop residues

We tested a simple aggregate indicator based on the stage of maturity to classify crop residues with respect to their potential to induce N 2 O emissions. The maturity stage, used to categorize crop residues as mature and immature residues, showed a large impact on soil N 2 O emissions. Immature residues, with an overall composition of high content of WSC, soluble NDS fraction and N returned, and low cellulose, hemicellulose, lignin, LCI and C: N ratio, stimulated the effect of residues on N 2 O emissions. Immature residues were composed by cover crops, vegetable residues and grassland renewal. Incorporation of mature residues on the other hand, had a weak influence on N 2 O emissions. Mature residues included a wide range of cereal residues in the form of straw. This simple and yet robust categorization offers a tool for predicting N 2 O emissions from crop residues obtained under a broad range of pedoclimatic conditions, encompassing intrinsically the different crop properties and management practices associated with their cultivation. Previous syntheses with a similar focus advocated for the use of residue C:N ratio as the best predictor of N 2 O emissions from crop residues [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF][START_REF] Muhammad | Regulation of soil CO2 and N2O emissions by cover crops: a meta-analysis[END_REF][START_REF] Shan | Effects of crop residue returning on nitrous oxide emissions in agricultural soils[END_REF]). Yet, our random-meta-forest approach revealed that the C:N ratio was not among the most solid predictors of residuederived N 2 O emissions (Supp. Fig. 1). This may be because C:N ratio only accounts for the balance between residue N mineralization/immobilization, whereas maturity class, by integrating a wide range of biochemical properties, also addresses other controlling factors of N 2 O emissions such as C availability and soil O 2 consumption. Therefore, the use of our simple classification according to maturity class provides new opportunities to identify crop residues requiring targeted agronomic practices for reducing N 2 O emissions, inform agricultural policy measures designed to mitigate climate change, and predict the potential impacts of crop rotations on N 2 O emissions.

Impact of abiotic conditions

Our random-meta-forest approach indicated that soil texture, soil bulk density, mean annual precipitation, mean annual temperature, and study length had little predictive power regarding N 2 O emissions from crop residues. Our hierarchical approach, using only data from field experiments that are representative for on-farm conditions, indicates that although some of these factors exert an influence on crop residue decomposition and associated N 2 O emissions (e.g. mean annual precipitation and soil texture; [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF], their overall importance is much lower than that of residue biochemical properties, and also lower than that of other abiotic conditions.

The most relevant soil and environmental factors among the ones tested in our study were soil pH, SOC, soil N and Aridity Index, with neutral soil pH, low SOC and low soil N content leading to lower N 2 O emissions. For the soil properties, these results are broadly in agreement with previous meta-analyses [START_REF] Chen | Soil nitrous oxide emissions following crop residue addition: a meta-analysis[END_REF][START_REF] Xia | Trade-offs between soil carbon sequestration and reactive nitrogen losses under straw return in global agroecosystems[END_REF][START_REF] Hu | The responses of soil N2O emissions to residue returning systems: a meta-analysis[END_REF]. The impact of SOC and soil N reflects their role as C and N sources for N 2 Oproducing microorganisms; low soil pH (acidic) increases N 2 O fluxes due to impairment of periplasmatic N 2 O reductase [START_REF] Bergaust | Denitrification response patterns during the transition to anoxic respiration and posttranscriptional effects of suboptimal ph on nitrogen oxide reductase in paracoccus denitrificans[END_REF][START_REF] Bakken | Regulation of denitrification at the cellular level: a clue to the understanding of N2O emissions from soils[END_REF]. While the importance of soil parameters has been summarized in previous meta-analyses, our study is the first one showing the role of the Aridity Index for crop residue-derived N 2 O emissions. The effect size of applying crop residues in humid climates was 2 times greater than in dry climates. As posited by [START_REF] Rochette | Estimation of N2O emissions from agricultural soils in Canada. I. Development of a country-specific methodology[END_REF], using annual precipitation can lead to biased conclusions because sites with precipitation exceeding 1000 mm are usually only found under humid climates, which could be tropical or temperate. The ratio between precipitation and evapotranspiration corrects for this bias, which explains why the Aridity Index was a better predictor of N 2 O emissions from crop residues than annual precipitation and annual temperature alone, and therefore its use should be promoted.

Implications

The strong positive effect of immature crop residues on N 2 O emissions has important implications for current efforts to achieve sustainable cropping systems, as residues falling in this category originate from cover crops, grasslands, and vegetables. Cover crops provide a wide range of ecosystem services [START_REF] Haruna | Cover crop management effects on soil physical and biological properties[END_REF], including reductions of N losses in the form of nitrate leaching [START_REF] Abdalla | A critical review of the impacts of cover crops on nitrogen leaching, net greenhouse gas balance and crop productivity[END_REF]. However, the increase in N 2 O emissions after their incorporation may partly compromise these benefits. Our analysis only included cover crop effects on N 2 O emissions after incorporation, but it is possible that cover crops may reduce soil N 2 O emissions during their growing phase by reducing soil mineral N availability, thereby countering their afterlife effects [START_REF] Han | N2O emissions from grain cropping systems: a meta-analysis of the impacts of fertilizer-based and ecologically-based nutrient management strategies[END_REF]. Including temporary grasslands in rotation with annual crops may reduce water pollution [START_REF] Parish | Multimetric spatial optimization of switchgrass plantings across a watershed[END_REF], soil erosion [START_REF] Pimentel | Annual vs. Perennial grain production[END_REF], and pests [START_REF] Werling | Perennial grasslands enhance biodiversity and multiple ecosystem services in bioenergy landscapes[END_REF], and promote biodiversity [START_REF] Meehan | Bird communities in future bioenergy landscapes of the upper Midwest[END_REF] and soil C sequestration [START_REF] Beniston | Soil organic carbon dynamics 75 years after land-use change in perennial grassland and annual wheat agricultural systems[END_REF]). Yet, our results suggest that the increased N 2 O emissions after grassland termination may largely offset the GHG balance benefits of increased soil C storage. Regarding crop residues of vegetables, it is worth noting that global vegetable production is rapidly expanding due to health guidelines advocating for an increase in vegetable consumption [START_REF] Norris | Alternative management practices improve soil health indices in intensive vegetable cropping systems: a review[END_REF]. However, the high N 2 O emissions induced by their residues (and during the vegetation period; [START_REF] Qasim | Global greenhouse vegetable production systems are hotspots of soil N2O emissions and nitrogen leaching: a meta-analysis[END_REF], and the high requirements for fertilizers, irrigation, and tillage needed for vegetable production call for an increased research effort to improve the environmental sustainability of these systems.

Our findings also have major implications for calculating N 2 O emission factors from crop residues. First, the relationship between supplied N with crop residues and N 2 O emissions is not linear (Fig. 2), which challenges the use of an emission factor that assumes linearity between N returned in residues and emitted N 2 O as assumed by the IPCC Tier 1 approach (IPCC, 2019). Second, using one single factor as residue N is not appropriate to explain the magnitude of the response of N 2 O emissions to the application of crop residues across multiple sites, years and cropping systems. The use of integrative quality criteria such as maturity class or residue type, are more promising approaches to constrain the variability regarding N 2 O emissions from crop residues. This classification could therefore be used to improve the estimation of N 2 O emission factors from crop residues, and it may be easily incorporated in the IPCC guidelines to estimate N 2 O emissions from agricultural sources.

CRediT authorship contribution statement

Diego Abalos: conceptualization, formal analysis, writing -review and editing; Tatiana Rittl: data curation, writing -original draft, review and editing; Sylvie Recous: conceptualization, data curation, writing -review and editing; Pascal Thiébeau: investigation, resources, data curation; Cairistiona F. E. Topp: data curation, review and editing; Kees Jan van Groenigen: formal analysis, review and editing; Klaus Butterbach-Bahl: review and editing; Rachel E. Thorman: data curation, review and editing; Kate E. Smith: review and editing; Ishita Ahuja: data curation, review and editing; Jørgen E. Olesen: project management, review and editing; Marina A. Bleken: review and editing; Robert M. Rees: review and editing; Sissel Hansen: conceptualization, coordination, writing -review and editing.

Fig. 1 .

 1 Fig. 1. Relative increase in field N 2 O emission (mean ± 95% CIs) with crop residue retention, as affected by crop type. The number of observations and studies are shown in parentheses. See section 2.1.2. for detailed descriptions of crop categories.

Fig. 2 .

 2 Fig. 2. Meta-analytic scatterplots of important factors controlling the effect of crop residue retention on N 2 O emissions (identified by the random-meta-forest approach); a) N returned with crop residues (kg N ha -1 yr -1 ), b) residue hemicellulose concentration (% DM), and c) residue cellulose concentration (% DM).

Fig. 3 .

 3 Fig. 3. Meta-analytic scatterplots of important factors controlling the effect of crop residue retention on N 2 O emissions (identified by the random-meta-forest approach) and influenced by N fertilizer addition: Water soluble carbon (% total C), soluble NDS fraction (% DM), and lignin (% DM). Graphs in the top row (blue; a, c, e) show results for studies without N fertilizer application; the bottom row (yellow; b, d, f) shows results for studies with N fertilizer application. Note that the scale of the horizontal axes for lignin are different for fertilized and unfertilized meta-regressions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Relative increase in field N 2 O emission (mean ± 95% CIs) with crop residue incorporation, as affected by residue maturity class. The number of observations and studies are shown in parentheses. See section 2.1.3. for a detailed description of the maturity criteria.

Fig. 5 .

 5 Fig. 5. Principal component analysis of the plant biochemical characteristics according to a) maturity class, and b) crop type. Ellipses indicate 95% confidence intervals.

Table 1

 1 Overview of factors used to categorize the crop residues included in our analysis. crops; (ii) pre-harvest destruction of the green haulm of potatoes, carrots, sugar beet and other root crops; (iii) harvest of brassicas and other vegetable crops (e.g., cauliflowers, cabbages and green salads), where leafy material is left in the field; (iv) ploughing down of herbage of grass leys

	Groups	Factors	Number of	Study categories and range for continuous factors
			observations
	Residue-N	N returned with crop residues (kg N ha -1 yr -1 )	297 a	5 to 418
	Crop type	Crop type	367	Cereal, cover crop, grassland, grain legumes, rice, sugar cane, vegetable, double cropping
	Residue type	Type of residue generated	367	Green plant biomass, mature above ground biomass, senescent plant biomass, straw
		Maturity index	367	Immature, mature
	Residue quality	Residue C:N ratio	357 a	<20, 20-60, >60
		Soluble NDS (% total DM)	344 a	6 to 71
		Cellulose (% total DM)	346 a	10 to 49
		Hemicellulose (% total DM)	345 a	6 to 55
		Lignin (% total DM)	352 a	2 to 26
		Lignocellulose Index (LCI)	344 a	0.041 to 0.306
		Water soluble carbon (% total C)	288 a	3 to 68
	Soil properties	Soil texture	222	Clay, loam, sandy
		Clay (%)	219	3.1 to 71%
		Soil pH	288 a	Acidic (<6), neutral (6-7), alkaline (>7)
		Soil organic carbon (g C kg -1 SDW)	310 a	2 to 98
		Soil bulk density (g cm -3 )	178	0.76 to 1.7
		Soil total N (g N kg -1 SDW)	262	0.14 to 6.4
	Climatic conditions	Annual mean precipitation (mm)	367 a	350 to 2115
		Annual mean temperature (°C)	227	5.3 to 27.4
		Aridity index	361 a	<1, >1
		Normalized days at 15 °C	192	4.7-1001

a Factors included as potential predictors in the random-meta-forest approach. catch

Acknowledgments

This study was supported by the FACCE ERA-GAS Cofund (European Research Areafor mitigation of Greenhouse gases from Agri-and Silviculture; id number: 696356), project ResidueGas (https://projects.au.dk/ residuegas/). The Danish funding was received from Innovation Fund Denmark (project 7108-00005B). The French funding was received from the Agence Nationale de la Recherche (ANR-17-EGAS-0003) and the National Institute for Food, Agriculture and the Environment (INRAE). The authors acknowledge funding in the UK from the Department for Environment, Food and Rural Affairs (Defra) and the Scottish Government's Strategic Research Programme. The authors acknowledge funding in Norway from the Research Council of Norway (Project No.: 276389). We would like to thank the researchers who were contacted and kindly supplied any missing information: Liz Baggs, Søren O. Petersen, Claudia Wagner-Riddle, Yuefang Zhang, Pierre-André Jacinthe, Caroline Buchen, Kari Dunfield, Maria Ernfors, Katharina Kesenheimer, Xiaoxi Li, Marco Siqueira Neto, Thorsten Reinsch, Reiner Ruser, Reinhard Well, Khagendra Baral, and Jaber Rahimi. The Chamber of Agriculture Lower-Saxony contributed with climatic data. Special thanks to Christian Werner for extraction of geodata.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi. org/10.1016/j.scitotenv.2021.152532.