Baoling Guo
email: baoling.guo@g2elab.grenobble-inp.fr

Seddik Bacha

Mazen Alamir

Amgad Mohamed

Cédric Boudinet

LADRC applied to variable speed micro-hydro plants: experimental validation

Keywords: Micro-hydro plants, variable speed, disturbances, LADRC, torque observer, robustness

To track the maximum power of a hydraulic turbine, the rotational speed needs to be adapted to flow rates variations, enabling higher speed tracking demands. Moreover, a Variable Speed Micro-Hydro Plant (VS-MHP) is a complex nonlinear system affected by large uncertainties. Active Disturbance Rejection Control (ADRC) can estimate and compensate the total disturbances in real-time, achieving higher robustness. A Linear-ADRC (LADRC) with the torque compensation is implemented to the speed control loop of a VS-MHP, which is tested under a real-time power Hardware-In-the-Loop benchmark. Experimental results prove its advantages for VS-MHP applications.

Introduction

Hydropower is widely used in the renewable energy context [START_REF] Ansel | Modelling and simulation of an autonomous variable speed micro hydropower station[END_REF][START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF][START_REF] Márquez | Dynamic modeling, simulation and control design of an advanced micro-hydro power plant for distributed generation applications[END_REF][START_REF] Guo | Variable speed micro-hydro power generation system: Review and Experimental results *[END_REF][START_REF] Valavi | Variable-speed operation of hydropower plants: Past, present, and future[END_REF]. Limited to the capacity of the Power Electronics (PE), variable speed technique is mainly implemented to small-or micro-hydro plants, which are mostly based on the "Run-of-River" (RoR). The RoR scheme usually has a low head and high water flow rates, Kaplan or Semi-Kaplan turbines become the most attractive prime mover for their high efficiency under such conditions 5 [START_REF] Ansel | Modelling and simulation of an autonomous variable speed micro hydropower station[END_REF][START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF]. The hydraulic efficiency of Kaplan-type turbine is sensitive to water flow rate variations and rotational speed changes, the variable speed operation can help find the best efficiency point when the hydraulic conditions change [START_REF] Guo | Variable speed micro-hydro power generation system: Review and Experimental results *[END_REF]. Besides, variable speed techniques can optimise the operation process by mitigating the cavitation effects [START_REF] Valavi | Variable-speed operation of hydropower plants: Past, present, and future[END_REF], alleviating water hammer disturbances [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF], and optimizing transient processes [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF]. Thus, Variable Speed Micro-Hydro Plants (VS-MHP) equipped with PE units gives rise to advantages in the hydraulic context [START_REF] Ansel | Modelling and simulation of an autonomous variable speed micro hydropower station[END_REF][START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF][START_REF] Márquez | Dynamic modeling, simulation and control design of an advanced micro-hydro power plant for distributed generation applications[END_REF][START_REF] Guo | Variable speed micro-hydro power generation system: Review and Experimental results *[END_REF][START_REF] Valavi | Variable-speed operation of hydropower plants: Past, present, and future[END_REF][START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF]. In 10 this paper, a VS-MHP of semi-Kaplan turbine coupled to a Permanent Magnet Synchronous Generator (PMSG) is considered. The hydropower generation system is connected to the power grids via back-to-back full PE converters.

The rotational speed of hydraulic turbine needs to be adapted to the water flow rates variations in order to achieve the Maximum Power Point Tracking (MPPT) in a VS-MHP, which enables higher speed tracking demands [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF]. MPPT operation provides an optimal speed reference for the speed control loop. Then, Proportion-Integration 15 (PI) controllers are commonly used for the speed regulation due to their simplified structure, high steady-state accuracy, and stability [START_REF] Jung | Adaptive PID Speed Control Design for Permanent Magnet Synchronous Motor Drives[END_REF]. Nevertheless, a VS-MHP is a nonlinear system disturbed by various internal and external uncertainties [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF]. Several internal disturbances must be considered: A PMSG is faced with various uncertainties such as parameters variations caused by the aging factor, the non-constant friction, and the nonlinear magnetic field effect [START_REF] Li | Adaptive speed control for permanent-magnet synchronous motor system with variations of load inertia[END_REF]; Additionally, PE itself is full of uncertainties including non-constant parameters caused by the aging factor, nonlinear variations due to heat effects. External disturbances mainly involve the following aspects:

The water dynamics inside the penstocks and the actuator's dynamics cause pressure head oscillations, which would result in torque fluctuations [START_REF] Mesnage | Modélisation et contrôle avancé pour les centrales de turbinage de moyenne et haute chute[END_REF]. Also, there are cases where the turbine operates in unstable regions, being unavoidable and resulting in severe vibrations [START_REF] Mesnage | Hydraulic-turbine start-up with "s-shaped" characteristic[END_REF]. The nonlinear features and large unpredictable disturbances in VS-MHPs would negatively affect the speed control performance. The increasing demands for high robustness and unavoidable limitations of traditional PI controllers have driven the engineers to apply more efficient control solutions [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF][START_REF] Jung | Adaptive PID Speed Control Design for Permanent Magnet Synchronous Motor Drives[END_REF][START_REF] Li | Adaptive speed control for permanent-magnet synchronous motor system with variations of load inertia[END_REF][START_REF] Kung | FPGA-Based Speed Control IC for PMSM Drive With Adaptive Fuzzy Control[END_REF][START_REF] Elghali | High-order sliding mode control of a marine current turbine driven doubly-fed induction generator[END_REF][START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF][START_REF] Guo | A review on ADRC based PMSM control designs[END_REF][START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF].

Large advanced control techniques such as adaptive Proportion-Differentiation-Integration (PID) [START_REF] Jung | Adaptive PID Speed Control Design for Permanent Magnet Synchronous Motor Drives[END_REF], Fuzzy control [START_REF] Kung | FPGA-Based Speed Control IC for PMSM Drive With Adaptive Fuzzy Control[END_REF], sliding mode control [START_REF] Elghali | High-order sliding mode control of a marine current turbine driven doubly-fed induction generator[END_REF], and Active Disturbance Rejection Control (ADRC) [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF][START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF][START_REF] Guo | A review on ADRC based PMSM control designs[END_REF][START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF] are proposed to improve the speed control performance. Due to high robustness and advantages of less depending on the knowledge of a plant model, ADRC is increasingly implemented to the speed control of electrical drive systems [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF][START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF][START_REF] Guo | A review on ADRC based PMSM control designs[END_REF][START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF]. The key module in the frame of ADRC is the Extended State Observer (ESO). The internal dynamics and external disturbances are lumped together as the generalized disturbances, which are defined as an extended state variable. Then, total disturbances are estimated via the ESO and actively compensated into the closed-loop dynamics in real-time, which achieves high robustness in dealing with uncertainties. Further, Linear-ADRC (LADRC) has been proposed in [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF], which make the controller design and parameter tuning more practical.

Large applications prove that it can achieve nearly the same performance as nonlinear ADRC with less complexity [START_REF] Li | Adaptive speed control for permanent-magnet synchronous motor system with variations of load inertia[END_REF][START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF][START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF][START_REF] Sun | A DSP-based active disturbance rejection control design for a 1-kW H-bridge DC-DC power converter[END_REF][START_REF] Liu | Active disturbance rejection controller based on permanent magnetic synchronous motor model identification and compensation[END_REF]. Moreover, the large hydraulic torque variations bring great estimation burden for ESO [START_REF] Guo | An anti-disturbance ADRC based MPPT for variable speed microhydropower plant[END_REF][START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF], which would affect the speed tracking performance. If the torque information could be obtained and appropriately used, the estimation efficiency would be enhanced. However, a physical torque sensor is still expensive; some sensorless approaches can observe the load torque precisely at the price of the extra computational complexity [START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF][START_REF] Zedong | A Rotor Speed and Load Torque Observer for PMSM Based on Extended Kalman Filter[END_REF].

A simple on-line torque identification approach discussed in [START_REF] Liu | Active disturbance rejection controller based on permanent magnetic synchronous motor model identification and compensation[END_REF] is employed in this design, which can track the hydraulic torque dynamics efficiently.

In this paper, a LADRC-based speed controller is designed to improve the speed tracking performance, thus achieving efficient MPPT operation. The key features of the proposed speed control strategy are as follow: 1) It achieves faster control dynamics and higher robustness against various external disturbances and internal uncertainties; 2) It introduces an on-line torque observer into the speed control loop, which enables more efficient ESO; 3) It is simple to design and tune the speed controller for practical engineers. This paper is organised as follows: In Section 2, the VS-MHP under study is modelled and described. In Section 3, a LADRC-based speed controller is designed, and a hydraulic torque observer is introduced to enhance the efficiency of ESO. In Section 4, the suggested LADRC-based speed controller is validated under a real-time power Hardware-In-the-Loop (HIL) benchmark. Conclusions are indicated in Section 5.

Variable speed micro-hydro plant

The VS-MHP under study includes three parts: Hydraulic turbine, PMSG, and Power Electronics interfaces, which are modelled and described in the following subsections. Besides, the employed MPPT technique and the control design of VS-MHP are explained before starting the LADRC-based speed controller design. In this paper, a propeller-type semi-Kaplan turbine of fixed blades and fixed guide vanes is considered, which has a similar operating feature to wind turbines [START_REF] Iman-Eini | Evaluation of loss effect on optimum operation of variable speed micro-hydropower energy conversion systems[END_REF]. The available hydraulic power P h (W ) depends on the net water head H(m) and the water flow rate Q (m 3 /s):

P h = ρ • g • H • Q (1)
where ρ (kg/m 3 ) is the density of water, g is the acceleration due to gravity (m/s 2 ).

Hydraulic turbine efficiency η strongly affects the net output mechanical power P m (W ) as follows:

P m = η • P h (2)
The hydraulic mechanical torque is represented by (3), and ω (rad/s) is the hydraulic turbine rotational speed.

T m = P m ω (3)
The steady-state efficiency of the semi-Kaplan hydraulic turbine is experimentally described in [START_REF] Márquez | Dynamic modeling, simulation and control design of an advanced micro-hydro power plant for distributed generation applications[END_REF], which is approximated by the following expressions:

η(λ, Q) = 1 2 [( 90
λ i + Q + 0.78)exp( -50 λ i )]3.33Q (4)
where λ i is given by

λ i = [ 1 λ + 0.089 -0.035] -1 (5)
where λ = RAω/Q, R (m) is the radius of the hydraulic turbine, A (m 2 ) is the area swept by rotor blades. The efficiency described by (4) and the generating hydraulic power curves (2) are shown in Figure 1. The figures 60 indicate that when water flow rates change, there exists an optimal rotational speed that maximizing the efficiency.

3

Mathematical model of PMSG

In the synchronous dq frame, the PMSG can be modelled as follows:

v d = L d di d dt + R s i d -ω e L q i q v q = L q di q dt + R s i q + ω e (L d i d + ψ f ) (6)
where v d and v q are the stator voltage vector of q-axis and d-axis, respectively; i d and i q are the stator current vector of q-axis and d-axis, respectively; L d and L q are the inductance of stator winding of q-axis and d-axis, respectively; R s is the stator winding resistance; ψ f is the permanent magnet flux vector, ω e is the electrical rotational speed.

65

Electromagnetic torque T e in the dq coordinate system is given by [START_REF] Jung | Adaptive PID Speed Control Design for Permanent Magnet Synchronous Motor Drives[END_REF], where n p is the number of pole pairs.

T e = 3 2 n p (ψ f i q + (L d -L q )i q i d ) (7)
If L d = L q or i d = 0, the expression ( 7) is simplified by

T e = 3 2 P n ψ f i q (8)
A mechanical dynamic model between hydraulic turbine and PMSG is expressed by

J dω dt = T e -T m -Bω (9)
where J is the total inertia; B is the friction factor; T m is the hydraulic driven torque.

Adaptive Perturb and Observe MPPT technique

Perturb and Observe (P&O) MPPT technique does not require the prior knowledge of the system or additional flow rate sensors. Such advantages make it preferable for VS-MHPs because it is not easy to measure flow rates exactly in real-time. The converter and generator losses can change the maximum power point location of the whole hydropower system [START_REF] Iman-Eini | Evaluation of loss effect on optimum operation of variable speed micro-hydropower energy conversion systems[END_REF]. Therefore, when the shaft speed is perturbed, the injected power to the grid is observed (but not the mechanical power of hydraulic turbine) [START_REF] Iman-Eini | Evaluation of loss effect on optimum operation of variable speed micro-hydropower energy conversion systems[END_REF]. The grid power and the speed variations are computed respectively by

∆P g = P g (k) -P g (k -1) ∆ω = ω(k) -ω(k -1) (10)
The perturbing direction for hydro turbine rotational speed is determined by

δ(t) = sgn(∆P g )sgn(∆ω). If δ(t)
is positive, the speed increases in order to make the system track towards the maximum power point; otherwise, the speed decreases [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF]. However, the speed tracking dynamics and oscillations appearing in the steady state highly depend on the perturbing coefficient K and the sampling time T e . To address the trade-off issue between steady oscillations and its tracking speed, an adaptive perturbing coefficient K adap is thus employed in [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF]. The speed reference is updated by

ω * (k) = ω(k) + t k-1 +Te t k-1 K adap δ(t)dt (11)
Finally, the output of the MPPT provides the optimal speed reference for the speed control loop.

Control of variable speed micro-hydro plant

The power architecture of VS-MHP with its control design is schematically presented in Figure 2. The PMSG is driven by a semi-Kaplan hydraulic turbine, and the full back-to-back PE converters are employed to transfer the hydraulic power as much as possible to the grid. The machine-side converter control consists of two cascaded control loops. Usually, the external speed control loop regulates the turbine shaft to run at an optimal rotational speed [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF]. As aforementioned, an adaptive MPPT technique is employed to generate the optimal speed reference. The rotational speed is determined by the dynamic function of the hydraulic torque and the machine electromagnetic torque. The machine torque is controlled through the stator current vector control upon a rotating dq synchronous reference frame. PI controllers are applied into the d-axis and q-axis current control loops, respectively. Further, a LADRC-based speed control with torque compensation is employed, it offers a fast dynamic response and high robustness under various disturbed conditions, which will be discussed in Section 3.

The grid-side converter control includes double closed control loops. The outer DC-link voltage loop ensures the stability of DC bus voltage. The inner loop needs to control the current to synchronize with the grid, achieving high quality injected current.

LADRC-based speed control strategy

The dual-loop control design of machine-side converter is presented in Figure 3. A first-order LADRC-based speed controller is introduced to improve the speed tracking performance, thus achieving efficient MPPT operation.

ESO can estimate external disturbances and internal dynamics in real-time. A proportional gain is employed to regulate the rotational speed. Note that the P&O MPPT operation imposes a speed tracking process before it reaches the maximum power point, which behaves like a transient process. The overshoot caused by the set-point input is mitigated, so the transient profile is neglected in this design. Before applying the LADRC paradigm, a key procedure is to reformulate a practical control issue to the canonical mode of cascade integrators as the expression [START_REF] Elghali | High-order sliding mode control of a marine current turbine driven doubly-fed induction generator[END_REF]. A systematic design procedure for a class of ADRC-based control applications is detailed in [START_REF] Guo | A review on ADRC based PMSM control designs[END_REF][START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF][START_REF] Han | From PID to active disturbance rejection control[END_REF]], which can be followed for the rotational speed control design.

                                         ẋ1 = x 2 . . . ẋ2 = x 3 ẋn-1 = x n ẋn = bu + f 0 (x 1 , x 2 , • • • , x n ) + f 1 (x 1 , x 2 , • • • , x n , w(t), t) ẋn+1 = h y = x 1 (12)
where x 1 , x 2 , • • • , x n are the state variables, y is the output, u is the control input, w(t) is the uncertain external disturbance, b is the control gain,

f 0 (x 1 , x 2 , • • • , x n ) represents the known disturbances, x n+1 = f 1 (x 1 , x 2 • • • , x n , w(t), t)
represents the unknown disturbances, which is defined as an extended stated, and

h = ḟ1 (x 1 , x 2 • • • , x n , w(t), t) is
the time derivative of disturbances [START_REF] Han | From PID to active disturbance rejection control[END_REF].

95

By simply manipulating expressions ( 8) and ( 9), the dynamic model can be rewritten by

dω dt = 3 2J n p ψ f i * q - T m J - Bω J (13)
Firstly, it is essential to choose the output y to be controlled and the input u to be manipulated. Based on [START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF],

ω is defined as the output y, i * q is chosen as the input u. Then, the order of LADRC-based controller is chosen according to the order of the plant model. A first-order LADRC is hence considered, and a second-order ESO is needed.

Further, a crucial step is how to reformulate the problem by lumping various known, unknown components into the total disturbances. The efficiency of LADRC largely depends on its key element ESO, and the accuracy of ESO is affected by the total disturbances [START_REF] Han | From PID to active disturbance rejection control[END_REF]. Expression [START_REF] Li | Output predictor-based active disturbance rejection control for a wind energy conversion system with pmsg[END_REF] indicates that speed variations and hydraulic torque fluctuations are the main external disturbances for the speed control loop. To mitigate the effects of large speed changes, -B ω/J is easy to be computed as the known disturbances, where ω can be estimated by ESO. In fact, the hydraulic torque fluctuations are the main external disturbances. The hydraulic torque information is needed to mitigate the effects of torque variations. While a physical torque sensor is expensive, an on-line torque identification approach based on the mechanical dynamic equation ( 9) is thus employed, which is detailed in the following part.

Now, let y = x 1 = ω, u = i * q , b 0 = 3 2J n p ψ f , x 2 = d total
is defined as an extended state variable, the time derivative of disturbances is given by ẋ2 = h. The speed control problem is formulated to a first-order standard model by [START_REF] Guo | A review on ADRC based PMSM control designs[END_REF], having the same form to [START_REF] Elghali | High-order sliding mode control of a marine current turbine driven doubly-fed induction generator[END_REF].

           ẋ1 = b 0 u + +f 0 (ω, Tm ) + d total ẋ2 = h y = x 1 (14)
where d total can be represented by [START_REF] Tong | Research on PMSM active disturbance rejection controller based on model compensation[END_REF]. It implies that, if Tm is well observed by the torque observer, the total disturbances estimated by ESO are very small.

d total (t) = (b -b 0 )i * q -( T m J + Bω J ) + ( Tm J + B ω J ) + w(t) (15)

Extended State Observer

Based on the canonical model ( 14), a second-order ESO is constructed as follows:

   ż1 = -β 1 (z 1 -ω) + b 0 i * q + z 2 + f 0 (ω, Tm ) ż2 = -β 2 (z 1 -ω) (16)
where z 1 , z 2 are the estimations of the output ω and the unknown disturbances d total , respectively; the positive gains β 1 = 2ω 0 and β 2 = ω 2 0 are chosen such that characteristic polynomial ( 17) is Hurwitz stable , and ω 0 represents the bandwidth of ESO [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF].

λ(s) = s 2 + β 1 s + β 1 2 = (s + ω 0 ) 2 (17)

Control law design

The estimated disturbances z 2 and the known disturbances f 0 (ω, Tm ) are both compensated in the feedback control channel, the feedback control law is represented by

i * q = k p (ω * -z 1 ) -[z 2 + f 0 (ω, Tm )] b 0 ( 18
)
where k p is the proportional gain of controller.

Within the ESO's bandwidth, the unknown disturbance is well estimated (z 2 ≈ d total ), the original disturbed plant model G p (s) is modified to be an integrator as follows [START_REF] Tatsumi | On the enhanced ADRC design with a low observer bandwidth[END_REF]:

Ĝp (s) = y(s) u 0 (s) ≈ 1 s (19)
Thus, a zero steady-state error for the speed can be achieved by employing a proportional controller. The closed-loop transfer function of the speed control loop is approximately represented by [START_REF] Iman-Eini | Evaluation of loss effect on optimum operation of variable speed micro-hydropower energy conversion systems[END_REF] where the closed-loop bandwidth of controller ω c ≈ k p .

H c (s) ≈ k p 1 s 1 + k p 1 s = 1 1 + s kp

Parameter tuning

Being different from a PI-type controller, LADRC-based controller is a Two-Degree-Of-Freedom (2DOF) controller, ω c and ω 0 can be adjusted independently [START_REF] Miklosovic | A robust two-degree-of-freedom control design technique and its practical application[END_REF]. Theoretically, the bigger the bandwidth, the faster the disturbance is observed and rejected by the controller [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF]. However, in a real application into a VS-MHP, the bandwidth is limited by sensor noises and practical considerations. An optimal bandwidth can be tuned under a real physical experimental benchmark, where the performance is maximized subject to its physical constraints [START_REF] Sun | A DSP-based active disturbance rejection control design for a 1-kW H-bridge DC-DC power converter[END_REF].

The tuning process is as follows:

-Firstly, ω c and ω 0 are designed to ensure the stability of the control system, and tested under the software environment. Generally, ω c is chosen to ensure the fast enough transient response, such as the setting time requirements [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF]. The bandwidth ω 0 is normally set to be 5-10 times of the controller's bandwidth ω c ; -Then, the control parameters are retuned and validated under the real-time physical experiments. The tests start with small values that can meet the basic performance requirements; -Further, both of them are incrementally increased until the control signal exceeds the noise tolerance of the q-axis current reference; -Finally, ω c and ω 0 are increased or decreased individually to find a best speed tracking performance.

In addition, PI-controller with the same bandwidth is designed to provide an equal comparison with the LADRCbased speed control. Details about adaptive MPPT design and PI-controller design can refer the work in [START_REF] Belhadji | Adaptive MPPT Applied to Variable-Speed Microhydropower Plant[END_REF][START_REF] Belhadji | Modeling and control of variable-speed micro-hydropower plant based on Axialflow turbine and permanent magnet synchronous generator (MHPP-PMSG)[END_REF].

Well-tuned control parameters are given in Table 1.

Hydraulic torque observer

The previous discussion in Section 3 suggests that the hydraulic torque information is essentially needed for an efficient LADRC-based speed controller, thus an on-line torque observer is employed in this design. By taking use of the information of i q and ω, the torque can be identified based on the following expression:

T m = K e i q -Bω -J ω (21)
To attenuate the negative effects from measurements noises, a low pass filter 1 T0s+1 is applied to the output side, then Laplace transformation is conducted to obtain the transfer function as follows:

Tm (s) = [K e i q (s) -Bω(s) -Jsω(s)] 1 T 0 s + 1 (22)
The expression ( 22) can be manipulated as follows:

Tm (s) = K e i q (s) -Bω(s) + J T0 ω(s) -Jsω(s) -J T0 ω(s) T 0 s + 1 = K e i q (s) -B ω(s) T 0 s + 1 - J T 0 ω(s) = ω 1 (s) - J T 0 ω(s) (23)
where

B = B -J T0 , ω 1 (s) = Keiq(s)-B ω(s) T0s+1
.

An inverse Laplace transformation is applied to get the time domain representation [START_REF] Belhadji | Modeling and control of variable-speed micro-hydropower plant based on Axialflow turbine and permanent magnet synchronous generator (MHPP-PMSG)[END_REF]. The complete observer design diagram is schematically illustrated in Figure 4. Finally, the control design of the machine-side converter is presented in Figure 5. This paper is dedicated to the machine-side speed control, so the control diagram of grid-side converter is not described here.

Tm = ω 1 - J T 0 ω (24) 0 1 T ' B 1 s 0 J T q i ˆm T  1   e K 1  0 1 T

Experimental validation

The experiments have been performed inside the real-time power HIL platform [START_REF] Gaztanaga | Real-Time Analysis of the Control Structure and Management Functions of a Hybrid Microgrid System[END_REF][START_REF] Munteanu | Hardware-in-the-loop-based simulator for a class of variable-speed wind energy conversion systems: Design and performance assessment[END_REF], which has been adapted for variable speed hydraulic tests as shown in Figure 6. Due to high costs and uncertain testing conditions of a real physical micro-hydro plant, a reduced-scale physical hydraulic simulator is developed in the laboratory. The hydro turbine model is implemented in the dSPACE hardware (DS1005), and its efficiency and water flow rates

Torque observer validation

Firstly, experiments are carried out to verify the effectiveness of the torque identification approach. Different types of torque disturbances are considered. In Figure 7 (a), the torque tracking performance of the starting process is presented. In Figure 7 (b), a step-up torque of 3N.m is first added and then it steps back to the original value.

T d (t) = D 3 H n (t) T11 (N (t), H n (t), γ(t)) (25a)
Q(t) = D 2 H n (t) Q11 (N (t), H n (t), γ(t)) (25b)
H n (t) = C h x h (t) (25c)
x h (t + τ ) = A h x h (t) + B h Q(t) (25d)
γ(t) = C a x a (t) (25e) x a (t + τ ) = A a x a (t) + B a γ c (t) (25f)
N (t + τ ) = N (t) + 30τ πJ t T d (t) (25g)
where the different parameters of the model given by equations (25a)-(25g) are defined in Table 2. 7 show that the torque observer is efficient to observe torque dynamics of the starting process, variations of step-up or step-down torque, and the cycle oscillated torque.

Variable speed control validation

To verify the speed tracking performance of LADRC-based controller, comparative experiments between strategies with torque compensation and without torque compensation have been conducted. To start the VS-MHP smoothly, a soft starting process is employed as shown in Figure 8. Firstly, a small initial torque is used to overcome the friction torque, the speed increases slowly; Then, the controller starts working under speed controlled mode in a no-load state; Finally, the MPPT mode is enabled and the rotational speed adapts to water flow rate variations. The results suggest that the overshoot of rotational speed is reduced with torque compensation, the total disturbances estimated z 2 is nearly zero, raising advantages when the torque is varying sharply.

Internal uncertainties test

The total inertia of VS-MHP can be expressed by J = J P M SG + J hdy , J P M SG is the inertia of the PMSG, and J hdy is the inertia of the hydraulic turbine. These two inertia values would vary when hydraulic conditions change [START_REF] Li | Adaptive speed control for permanent-magnet synchronous motor system with variations of load inertia[END_REF]. Suppose that the knowledge of system inertia is not precise, the controller is designed with different inertia values, but tested in the same benchmark. Figure 11 shows the speed tracking responses under different inertia values. It is expected that the increased inertia degrades the speed tracking performance, however, even the inertia value is increased to be 4 times of the normal value, the speed tracking process is stable and still under the toleration. The proposed LADRC-based speed control achieves high robustness in the presence of inertia variations.

Conclusions

In this paper, a VS-MHP of semi-Kaplan turbine coupled to a PMSG is studied. The hydropower generation system is connected to the power grids via back-to-back full PE converters. A VS-MHP is a classical nonlinear system disturbed by large uncertainties. ADRC-based controllers can estimate and compensate the influences of internal and external disturbances in real-time, which achieves high robustness and fast control dynamics. Thus, a LADRC-based speed controller is designed to improve the speed tracking performance, enabling efficient MPPT operation. To mitigate the negative effects caused by hydraulic torque variations, a torque observer is incorporated into the speed-control loop.

Different operation schemes are experimentally validated under a real-time power HIL benchmark. Some conclusions can be indicated: Firstly, the on-line torque observer can efficiently observe various torque disturbances such as torque dynamics of the starting process, step variations, and cycle oscillated torques. Further, the rotational speed can fastly track the speed reference during the starting process and the MPPT mode as well. The speed overshoot is reduced by compensating the hydraulic torque, particularly at the point when water flow rates vary.

Robustness tests are carried out considering external disturbances and internal dynamics as well. Compared with PI-controllers, the proposed LADRC-based speed controller achieves higher robustness and faster dynamics in the presence of various torque disturbances and system inertia variations. Further, the disturbance rejection ability of speed control is improved by introducing the torque observer. Experimental results prove its advantages for VS-MHP applications.

2 2. 1 .

 21 Mathematical model of hydraulic turbine

Figure 1 :

 1 Figure 1: Variable speed Semi-kaplan turbine model, with fixed head of H = 1m, R = 0.25m

Figure 2 :

 2 Figure2: Power architecture of VS-MHP and its global control design, where i (abc) represents the components of the machine current vector, Tm is the estimated hydraulic torque, V dc is the DC-link voltage, i g(abc) represents the components of the grid current vector, v g(abc) represents the components of the grid voltage vector, Lg and Rg are the inductance and resistance of the grid filter, respectively.

Figure 3 :

 3 Figure 3: Dual-loop control design of machine-side converter, Ke = 3 2 npψ f the torque constant, n the measurement noises

Figure 4 :

 4 Figure 4: Schematic of the hydraulic torque observer

  135

Figure 5 :

 5 Figure 5: Control design of machine-side converter, where θe is the electrical angle phase, v * a , v * b , v * c are the machine voltages references, Sa, S b , Sc are the switching signals, SVPWM: Space Vector Pulse Width Modulation.

Figure 6 :

 6 Figure 6: Real-time HIL-based VS-MHP experimental benchmark

Table 1 : 5 ,

 15 Parameters of the experiment benchmark Parameters of experimental benchmark Kmax = 2, Te = 0.1s kp = 2.5, ki = 333 ωc = 30, ω0 = 150

Figure 7 :

 7 Figure 7: Torque observer performance under different types of hydraulic torque, T * m is the emulated torque reference

Figure 8 :Figure 9

 89 Figure 8: Starting process of the VS-MHP, with ω the rotational speed (20rad/s/div), z 1 the estimation of the rotational speed (20rad/s/div), z 2 the estimation of total unknown disturbances (100/div)

Figure 9 :

 9 Figure 9: MPPT performance when water flow rates vary, with ω the rotational speed (20rad/s/div), z 1 the estimation of rotational speed (20rad/s/div), z 2 the estimation of total unknown disturbances (100/div)

4. 3 . 1 .Figure 10 :

 3110 Figure 10: Robustness of speed control under torque disturbances: (a) Step-up torque disturbance, (b) Step-down torque disturbance, (c) Oscillated torque disturbance, with ωerr = ω -ω *

Figure 11 :

 11 Figure 11: Performance comparisons under different inertia values

Table 2 :

 2 The definitions of the hydraulic dynamic variables Water flow rate at the side of the turbine m 3 /s A h , B h , C h Parameters of the penstocks model Q11 Regression model of the flow rate Hill charts A a , B a , C a Parameters of the actuator model Experimental results in Figure

	T d	Oscillated hydraulic torque (N.m)	x h	Internal states of the penstocks
	D	Turbine's diameter (m)	x a	Actuator's internal states
	H n	Net pressure head (m)	γ	Turbine's guide vanes opening (ratio)
	T11 Regression model of the torque Hill charts	γ c	Guide vanes angle set-point (ratio)
	N	Turbine's rotational speed (rpm)	J t	Turbine's moment of inertia (Kg.m 2 )
	Q

n p ψ f is the torque constant.

Acknowledgement

This paper is supported by PSPC Innov'hydro project, which specifically brings together GE (General Electric) Renewable, EDF (Electricity of France), Grenoble INP and other key players in the hydroelectric sector, in France.

Further, torque oscillations happen due to pressure head oscillations, a cycle oscillated torque is added to the basis torque as shown in Figure 7 (c). The numerical model of oscillated torque is described in the following part.

The torque oscillations result from the water dynamics inside the penstocks and the actuator's dynamics [START_REF] Mesnage | Modélisation et contrôle avancé pour les centrales de turbinage de moyenne et haute chute[END_REF].

The numerical model of torque oscillations is given by equations (25a)-(25g) [START_REF] Mesnage | Modélisation et contrôle avancé pour les centrales de turbinage de moyenne et haute chute[END_REF]. The mathematical models of the water dynamics inside the penstocks are given by equations (25c)-(25d). The dynamic actuator considered


Laure Soucémarianadin
email: lauresoucémarianadin@geologie.ens.fr

Moritz Reisser

Lauric Cécillon

Pierre Barré

Manuel Nicolas

Samuel Abiven

Pyrogenic carbon content and dynamics in top and subsoil of French forests

Keywords: pyrogenic carbon, benzene polycarboxylic acids, soil type, solubilization, subsoil, forest soils

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

Soil pyrogenic carbon (PyC) represents about 10-15% of total soil organic carbon (SOC) worldwide [START_REF] Forbes | Formation, transformation and transport of black carbon (charcoal) in terrestrial and aquatic ecosystems[END_REF][START_REF] Preston | Black (pyrogenic) carbon in boreal forests: A synthesis of current knowledge and uncertainties[END_REF][START_REF] Bird | The pyrogenic carbon cycle[END_REF]Reisser et al., 2016 and references therein) and can locally reach up to 35-45% [START_REF] Schmidt | Charred organic carbon in German chernozemic soils[END_REF][START_REF] Skjemstad | Charcoal carbon in U.S. agricultural soils[END_REF].

Deep (> 30 cm) SOC accounts for a significant amount of total SOC [START_REF] Jobbágy | The vertical distribution of soil organic carbon and its relation to climate and vegetation[END_REF], and is older and more stable than surface SOC [START_REF] Rumpel | Deep soil organic matter---a key but poorly understood component of terrestrial C cycle[END_REF][START_REF] Balesdent | Atmosphere-soil carbon transfer as a function of soil depth[END_REF]. This raises the question of how much PyC contributes to deep SOC stocks. However, to date, studies on PyC have mostly focused on surface soils, where the PyC is deposited after a fire. Only a few studies report PyC in deeper soil horizons (e.g., [START_REF] Rodionov | Black carbon in the zonal steppe soils of Russia[END_REF][START_REF] Guggenberger | Storage and mobility of black carbon in permafrost soils of the forest tundra ecotone in northern Siberia[END_REF][START_REF] Rodionov | Black carbon in grassland ecosystems of the world[END_REF][START_REF] Vasilyeva | Pyrogenic carbon quantity and quality unchanged after 55 years of organic matter depletion in a Chernozem[END_REF][START_REF] Abney | Post-wildfire erosion in mountainous terrain leads to rapid and major redistribution of soil organic carbon[END_REF] and out of the 55 studies considered in the meta-analysis of [START_REF] Reisser | Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks[END_REF], only 4 considered a complete soil profile.

Since PyC has to be transferred from the surface to reach deeper soil layers, its quality may be prone to changes with depth that have not been studied yet.

Drivers explaining PyC stocks in soil have been often related to fire properties. However, in a recent global study, [START_REF] Reisser | Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks[END_REF] showed that soil properties like pH or clay content had a stronger impact on the PyC concentration than fire characteristics, land use or climate. It is yet unclear if this global picture still holds at a landscape or regional scale.

In this study, we took advantage of a large network of forest sites and an established PyC characterization and quantification method to address the following research questions:

• Does PyC contribute significantly to total SOC in the subsoil?

• Does PyC quality evolve with depth?

We explored the effects of depth and of two potential drivers (soil class and vegetation type) on PyC quantity and quality using a balanced design.

We analysed soil samples collected from sites of the French national network for the long term monitoring of forest ecosystems (''RENECOFOR''). This network consists of 102 sites in managed forests, with even-aged stands. Mineral soil samples of the network were collected by layers up to 1-m depth to ensure comparability among sites and facilitate subsequent sampling to monitor SOC evolutions. Here, we chose to focus on samples collected at two depths: 0-10 cm, hereafter topsoil and 40-80 cm, hereafter subsoil. Details related to sample collection are available in Soucémarianadin et al. (2018a). We selected 22 sites, out of the entire network, to cover a diversity of soil classes and vegetation types. Soil classes were the following: (i) entic Podzols, (ii) eutric Cambisols and (iii) dystric Cambisols (IUSS Working Group, 2015). Forest vegetation was separated into two types with stands dominated by either (i) coniferous [silver fir (Abies alba Mill.) or Norway spruce (Picea abies (L.) H. Karst.)] or (ii) deciduous [beech (Fagus sylvatica L.)] trees. Study plots were mainly found in midslope or upslope positions (Table S1).

Bulk < 2 mm-sieved soil samples were ground and carbon concentrations were determined by dry combustion. Samples with carbonates were first decarbonated [START_REF] Harris | Acid fumigation of soils to remove carbonates prior to total organic carbon or CARBON-13 isotopic analysis[END_REF]. PyC was measured using the benzene polycarboxylic acids (BPCA) molecular marker method, which consists of extracting organic molecular markers specific from PyC [START_REF] Glaser | Black carbon in soils: The use of benzene carboxylic acids as specific markers[END_REF][START_REF] Brodowski | Revised black carbon assessment using benzene polycarboxylic acids[END_REF]. We followed the protocol described in [START_REF] Wiedemeier | Improved assessment of pyrogenic carbon quantity and quality in environmental samples by high-performance liquid chromatography[END_REF]. We used a conversion factor of 2.27 [START_REF] Glaser | Black carbon in soils: The use of benzene carboxylic acids as specific markers[END_REF] to calculate PyC proportion (% SOC) and content (g PyC•kg -1 soil). While this ratio has been criticized in the literature (e.g., [START_REF] Schneider | The benzene polycarboxylic acid (BPCA) pattern of wood pyrolyzed between 200 °C and 1000 °C[END_REF], we used it to compare our values with previous reports of PyC quantitative estimates [START_REF] Reisser | Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks[END_REF]. We also used the distribution patterns of BPCA markers as a quality indicator of the PyC. High proportion of B5CA and B6CA were indicative of more condensed and aromatic material [START_REF] Wiedemeier | Aromaticity and degree of aromatic condensation of char[END_REF], while B3CA and B4CA were the product of small condensed units of 3 aromatic rings minimum [START_REF] Ziolkowski | Quantification of black carbon in marine systems using the benzene polycarboxylic acid method: A mechanistic and yield study[END_REF].

We used paired t-test and multivariate models (see Soucémarianadin et al. (2018b) for details) to assess the effects of the three factors (depth, soil class and vegetation type) on the quantity and quality of PyC using R software [START_REF] Core | R: A language and environment for statistical computing[END_REF]. A generalized least squares function [START_REF] Pinheiro | Nlme: Linear and Nonlinear Mixed Effects Models[END_REF] was used for the linear mixed models, with forest site as a random effect. Model selection was implemented with a top-down strategy.

The mean PyC content in these French forest topsoils was 3.1 g PyC•kg -1 soil, representing 4.4% of total SOC (Table 1). This proportion was comparable to topsoils in Switzerland (0.6-4.7% of total SOC; Reisser, 2018) and in the "Centre" region of France (mean PyC = 5.3% of total SOC; [START_REF] Paroissien | Mapping black carbon content in topsoils of central France[END_REF]. However, our average PyC proportion to SOC was relatively low when compared to global values that average 13.7% of total SOC [START_REF] Reisser | Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks[END_REF]. This could be linked to the ongoing fire suppression over Europe that started in the 18 th century [START_REF] Pyne | Vestal Fire: An Environmental History, Told Through Fire, of Europe and Europe's Encounter with the World[END_REF]. The mean distribution of total BPCA among the four markers (Fig. 1; Table S2) was quite similar to what have been measured in Swiss forest topsoils [START_REF] Reisser | Fire's black legacy: The fate of pyrogenic carbon in soils across scales[END_REF]. In these topsoils, B5CA represented a proportion of total BPCA (23-51%; Reisser, 2018) similar to the one we observed in French forest soils. We also observed similar B6CA proportions to those reported in topsoils of alpine Leptosols (about 25%; [START_REF] Eckmeier | Black carbon contributes to organic matter in young soils in the Morteratsch proglacial area (Switzerland)[END_REF]. However, a few samples had very different proportions of B5CA and B6CA (Fig. S1). Specifically, the PyC in HET09 topsoil had a very low B5CA proportion (6%; Fig. S1), while PyC in HET21 and SP25 subsoils displayed relatively low B6CA (5% and 1%, respectively) and high B5CA proportions (52% and 55%, respectively; Fig. S1). Moreover, given the relative proportion of B6CA, around 30% (Fig. 1; Table S2), we assumed that soot contribution to total PyC had been minimal (e.g., [START_REF] Li | Organic matter source and degradation as revealed by molecular biomarkers in agricultural soils of Yuanyang terrace[END_REF] and that the PyC present in these forest soils originated mostly from biomass burning through forest fire.

Historical records showed that charcoal had been produced on-site in some of these forests (Table S1). However, no trace of kilns was observed on the sampling plots and although anthropogenic charcoal production might have constituted a source of PyC in some of our study plots, it is thus unlikely to be significant.

Although absolute PyC content decreased significantly with depth, PyC proportion tended to increase with depth (Fig. 2; Table 1). This increase was however not significant (Table 2). These results suggest a preferential transfer of PyC to the subsoil compared to other SOC moieties. Studies including samples from both topsoil and subsoil (> 30-cm depth) are not common, as shown by the meta-analysis of [START_REF] Reisser | Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks[END_REF], but some previous observations agreed with our results. In Longleaf pine forests of the southern USA, despite PyC content decreasing significantly with depth in most sites, PyC proportion (5-7% of SOC) was independent of depth [START_REF] Butnor | Vertical distribution and persistence of soil organic carbon in fire-adapted longleaf pine forests[END_REF]. A tendency for PyC proportion to increase with depth was observed in Russian Chernozems (from 7% to 10% SOC in 0-10 cm and 30-50 cm, respectively), in a mixed-grass savanna in Texas (5-9% BPCA-C at 0-10 cm vs. 7-13% BPCA-C at 10-20 cm) and in black soils of the Southern Alps [START_REF] Dai | Black carbon in a temperate mixed-grass savanna[END_REF][START_REF] Hammes | Centennial black carbon turnover observed in a Russian steppe soil[END_REF][START_REF] Eckmeier | Preservation of fire-derived carbon compounds and sorptive stabilisation promote the accumulation of organic matter in black soils of the Southern Alps[END_REF].

However, non-linear evolutions of PyC proportion with depth have been reported [START_REF] Rodionov | Black carbon in the zonal steppe soils of Russia[END_REF][START_REF] Hammes | Centennial black carbon turnover observed in a Russian steppe soil[END_REF][START_REF] Soucémarianadin | Pyrogenic carbon stocks and storage mechanisms in podzolic soils of fire-affected Quebec black spruce forests[END_REF] and distribution of PyC proportion in the complete soil profile cannot thus be inferred from our results. Downwards transfer of PyC by leaching, either through particulate or soluble forms, has been evidenced in various field conditions [START_REF] Hockaday | Direct molecular evidence for the degradation and mobility of black carbon in soils from ultrahighresolution mass spectral analysis of dissolved organic matter from a fire-impacted forest soil[END_REF][START_REF] Leifeld | Mobility of black carbon in drained peatland soils[END_REF][START_REF] Major | Fate of soil-applied black carbon: Downward migration, leaching and soil respiration[END_REF][START_REF] Santos | Thermal alteration of water extractable organic matter in climosequence soils from the Sierra Nevada, California[END_REF]. Although PyC quality did not change drastically between the top and deep soil layers (Fig. 1), we observed a significant decrease in B3CA and a concomitant significant increase in more condensed B5CA (Table S2). Abiven et al. (2011) observed that, over time, B5CA tended to preferentially solubilize, contrary to B3CA. These quality changes were also observed when comparing dissolved organic matter and particulate organic matter in rivers [START_REF] Wagner | In-stream sources and links between particulate and dissolved black carbon following a wildfire[END_REF]. Our results therefore suggest a significant contribution of soluble PyC to the deeper soil horizons. Meanwhile, we did not observe any significant change in the B4CA/B6CA or B5CA/B6CA ratios with depth (Table S2). As these ratios are mainly affected by PyC decomposition in soils, this result suggests that no major PyC degradation occurred in our soil profiles [START_REF] Singh | Transformation and stabilization of pyrogenic organic matter in a temperate forest field experiment[END_REF][START_REF] Guggenberger | Storage and mobility of black carbon in permafrost soils of the forest tundra ecotone in northern Siberia[END_REF].

While PyC solubilization seems to be consistent with the BPCA distribution patterns, alternative transfer processes may have occurred. For example, mesofauna could contribute to PyC transfer through bioturbation, but this was probably limited in our coniferous sites and in our acidic soils, which tend to have few earthworms. Particle movement in the soil macro-porosity after PyC fragmentation could have also played an important role in plots where charcoal pieces have been observed (Table S1). Erosion could impact horizontal transfer of PyC in sloping landscapes (e.g., [START_REF] Rumpel | Preferential erosion of black carbon on steep slopes with slash and burn agriculture[END_REF][START_REF] Abney | Post-wildfire erosion in mountainous terrain leads to rapid and major redistribution of soil organic carbon[END_REF]. Further burial of PyC in depositional landforms could affect its vertical distribution [START_REF] Abney | Pyrogenic carbon erosion: Implications for stock and persistence of pyrogenic carbon in soil[END_REF]. However this vertical transfer process was unlikely to dominate in our study plots given their geomorphologic position (Table S1) and the absence of colluvial soils in our data set.

Multivariate models showed that PyC quantity (content) and quality (relative proportion of B3CA) were influenced by soil class, while vegetation type was only marginally influential in B6CA distribution (Table 2). PyC proportion was higher in Podzols subsoils than in both Cambisols subsoils (Fig. 2). This could be due to an increased solubilization in Podzols caused by the podzolization process itself, during which transfer of organic matter and sesquioxides takes place [START_REF] Lundström | The podzolization process. A review[END_REF][START_REF] Buurman | Podzolisation and soil organic matter dynamics[END_REF]. High dissolved PyC concentrations were indeed measured in the E horizon of Podzols from sites that experienced severe burning around 100 years prior to sampling [START_REF] Santos | Impact of a historical fire event on pyrogenic carbon stocks and dissolved pyrogenic carbon in Spodosols in Northern Michigan[END_REF]. Similar to OC, PyC was shown to be retained and to accumulate in podzolic B horizons [START_REF] Soucémarianadin | Pyrogenic carbon stocks and storage mechanisms in podzolic soils of fire-affected Quebec black spruce forests[END_REF][START_REF] Santos | Impact of a historical fire event on pyrogenic carbon stocks and dissolved pyrogenic carbon in Spodosols in Northern Michigan[END_REF]. In these horizons, PyC was associated with the fine fraction and PyC content was correlated with sesquioxide content [START_REF] Soucémarianadin | Pyrogenic carbon stocks and storage mechanisms in podzolic soils of fire-affected Quebec black spruce forests[END_REF]. PyC accumulation was also observed in other soil types rich in Fe and Al oxides [START_REF] Major | Fate of soil-applied black carbon: Downward migration, leaching and soil respiration[END_REF][START_REF] Cusack | Mineralogical controls on soil black carbon preservation[END_REF][START_REF] Rumpel | Contrasting composition of free and mineral-bound organic matter in top-and subsoil horizons of Andosols[END_REF]. These previous observations suggest that dissolved PyC is likely to be stabilized in the subsoil of podzolic soils either through organo-mineral or organometallic complexes.

The marginal effect of vegetation type on PyC content, proportion and quality compared with soil class may be linked to the time frame that needs to be considered for PyC evolutions. In these managed forests, trees were planted and current species may not reflect species composition over the last centuries, which is the time frame related to PyC turnover. Conversely, soil class encompassed long temporal scale (millennium) matching PyC longer turnover and infrequent inputs into the soil, resulting in patterns over longer timescale.

A few conclusions could be drawn from the observations we made in the soil profiles of these twenty-two French forest sites. First, PyC transferred to and accumulated in the subsoil in all soil types. Second, PyC appeared to transfer downwards preferably through solubilization, specifically in podzolic soils. Third, the accumulation of PyC in the subsoil was soil class dependent, but not related to present vegetation. Yet, our investigation presents certain limitations and more studies on PyC in the complete soil profile should be carried on in diverse soil types to further improve our understanding of PyC transfer and accumulation in mineral soils.

Tables and figures captions

Table 1. Range (mean, minimum, maximum) of PyC proportion (percentage of total SOC) and content (g PyC/kg soil), ratio of B3CA/B5CA markers, SOC content (g OC/kg soil), pH in water and texture for topsoil and subsoil layers in the 22 study plots. * data from [START_REF] Ponette | RENECOFOR -Chimie Des Sols Dans Les 102 Peuplements Du Réseau : Campagne De Mesures 1993-95[END_REF].

Fig. 1 .

 1 Fig. 1. Mean (± standard deviation) relative distribution of total BPCA among the four

Fig. 2 .

 2 Fig. 2. Comparison of PyC proportion (% of total SOC) in topsoil (0-10 cm) and subsoil

Table 2 .

 2 Details of models and their significant terms selected to explain variations in PyC quantity and quality in the 22 soil profiles. All models used a gls function (see details in the Statistical analyses section). NS (non signifiant): p > 0.1; T (trend) p < 0.1;

	* p < 0.05; ** p < 0.01 and *** p < 0.001.

Table 1 .

 1 Range (mean, minimum, maximum) of PyC proportion (percentage of total SOC) and content (g PyC/kg soil), ratio of B3CA/B5CA markers, SOC content (g OC/kg soil), pH in water and texture for topsoil and subsoil layers in the 22 study plots. . max. mean min. max. mean min. max. mean min. max. mean min. max. mean min. max. mean min. max.

	Factor	classes	PyC proportion (% SOC)	PyC content (g BPCA-C/kg soil)	B3CA/B5CA ratio	SOC content (g/kg soil)	pHwater	clay content (%)*		sand content (%)*
	DEPTH	n mean mintopsoil 22 4,4 1,9 10,3	3,1	0,8 10,4 0,22 0,12 0,41 72,4 22,3 147,1	5,0 4,1 7,5	27	6	45	32	1	83
		subsoil 22	8,8 1,4 44,3	1,2	0,1	6,1 0,16 0,09 0,2 13,4 3,3 42,4	5,9 4,6 8,7	20	5	53	42	6	85
		All 44	6,6 1,4 44,3	2,2	0,1 10,4 0,19 0,09 0,41 42,9 3,3 147,1	5,5 4,1 8,7	23	5	53	37	1	85

Table 2 .

 2 Details of models and their significant terms selected to explain variations in PyC quantity and quality in the 22 soil profiles. All models used a gls function (see details in the Statistical analyses section). NS (non signifiant): p > 0.1; T (trend) p < 0.1; * p < 0.05; ** p < 0.01 and *** p < 0.001.

	Response variable Transformation	Predictors in final model and level of significance	n a
	PyC proportion	1/	depth NS + soil NS + depth × soil T	44
	PyC content	1/sqrt	depth*** + soil**	44
	B3CA	1/	depth** + soil T	44
	B5CA	1/	depth**	43
	B6CA	no	veg T + soil NS + veg × soil T	42
	B4CA	1/	no significant factor	44

a outliers were identified with Cook's distance.

Acknowledgments 10 This work was supported by the French Environment & Energy Management Agency (ADEME) [APR REACCTIF, piCaSo project], Campus France [PRESTIGE-2015-3-0008], the SOCUTE project funded by the "Emergences" program from ville de Paris, National Research Program 68 "Soil as a Resource" of the Swiss National Science Foundation and the University of Zurich Research Priority Programme on 'Global Change and Biodiversity'. We acknowledge three anonymous reviewers for their time and effort in commenting on the manuscript.


Sébastien Galtier
email: sebastien.galtier@u-psud.fr

Sergey V Nazarenko
email: sergey.nazarenko@inphyni.cnrs.fr

Éric Buchlin
email: eric.buchlin@ias.u-psud.fr

Simon Thalabard
email: simon.thalabard@impa.br

Nonlinear Diffusion Models for Gravitational Wave Turbulence

Keywords:

A fourth-order and a second-order nonlinear diffusion models in spectral space are proposed to describe gravitational wave turbulence in the approximation of strongly local interactions. We show analytically that the model equations satisfy the conservation of energy and wave action, and reproduce the power law solutions previously derived from the kinetic equations with a direct cascade of energy and an explosive inverse cascade of wave action. In the latter case, we show numerically by computing the second-order diffusion model that the non-stationary regime exhibits an anomalous scaling which is understood as a self-similar solution of the second kind with a front propagation following the law k f ∼ (t * -t) 3.296 , with t < t * . These results are relevant to better understand the dynamics of the primordial universe where potent sources of gravitational waves may produce space-time turbulence.

I. INTRODUCTION

The nonlinear nature of the Einstein equations of general relativity suggests that space-time can be turbulent. Such a turbulence has been studied in the context of spinning black holes [1][2][3][4] by using the gravity-fluid correspondence. It is shown numerically that such a system can display a nonlinear parametric instability with transfers reminiscent of an inverse cascade; the precise mechanism is, however, not totally understood. The possibility of having a turbulent cascade within the metric perturbations was already discussed in the past [5] but it is only recently that a rigorous theory has been proposed for the regime of gravitational wave (GW) turbulence [6]. The presence of small nonlinearities has been exploited to derive a four-wave kinetic equation which describes the spectral transfers of energy and wave action. In other words, the theory explains the nonlinear evolution of weak ripples on the Poincaré-Minkowski flat space-time metric. The theory is limited to a 2.5+1 diagonal metric tensor [7], which means that we only have a two-dimensional spatial dependence with no restriction on the metric components; the time variable is an additional variable (hence, the notation +1). This theory includes only one type (+) of GW (diagonal contribution); the × GW corresponds to a non-diagonal metric perturbation which is excluded by construction. Besides the kinetic equations, the main results obtained are the derivation of its power-law (constant flux) solutions and the demonstration that we have a direct energy cascade and an explosive inverse cascade of wave action (a property of finite-capacity turbulence systems) with a priori the possibility to excite Fourier modes from the injection wavenumber k I to k = 0 in a finite time. However, and as discussed in [6], such a transfer driven by GW turbulence stops at an intermediate scale where the turbulence regime becomes strong. Note that the change of regime does not preclude the possibility to extend such an inverse cascade to k = 0 in a finite time [8]. Generally speaking, we may say that space-time turbulence is likely to be a relevant regime for describing the very early universe soon after the Planck's time when the universe emerges from a quantum foam [9].

The main goal of the present paper is to study the properties of GW turbulence further and to compare it with other turbulent systems. The kinetic equation derived in [6] is, however, too complicated for detailed analytical or numerical studies. Therefore, our strategy is to derive nonlinear diffusion models (also called differential models) which correspond to strongly local interaction approximations of the kinetic equation. This type of reduction, first introduced by Leith in 1967 [10] to study three-dimensional Navier-Stokes turbulence, is quite common and plays an important role for achieving a qualitative and even quantitative understanding of turbulence in various physical situations, for both weak and strong turbulence [11][12][13]. For example, the anomalous scaling of the non-stationary spectrum, a property of finite-capacity turbulence systems first detected in the numerical simulations of some kinetic equations [14,15], has been easily studied in other systems using differential models [16][17][18][19][20][21]. In particular, it was established that the anomalous exponent is independent of the initial conditions but varies for different models of the same physical system.

In the present paper, we first recall the main properties of the kinetic equation of GW turbulence (section II) derived in [6], and then we propose two diffusion models for GW turbulence: a fourth-order and a second-order model which are introduced in sections III and IV, respectively. Their derivation is based, in particular, on the phenomenology of wave turbulence that was introduced in [6]. Numerical simulations of the second-order diffusion model are then performed to study the form of the front propagation during the inverse cascade of wave action. The results are presented in section V. Conclusions are discussed in section VI.

II. KINETIC EQUATION

Let us consider the metric tensor g µν = η µν + h µν , with η µν the Poincaré-Minkowski flat space-time metric and |h µν | ∼ h 1 a weak perturbation. Then, the threedimensional wave action spectrum n(k) scales like h 2 /k and the kinetic equation of GW turbulence writes [6]

∂n(k) ∂t = 4π |T kk3 k1k2 | 2 δ(ω k + ω k3 -ω k1 -ω k2 ) (n k1 n k2 n k3 + n k n k1 n k2 -n k n k2 n k3 -n k n k1 n k3 ) δ(k + k 3 -k 1 -k 2 ) dk 1 dk 2 dk 3 , (1)
where ω k = c|k| = ck (c is the speed of light) and T kk3 k1k2 is the coefficient describing the 4-wave interactions; it is characterized by the homogeneous property

T ξk ξk3 ξk1 ξk2 = T kk3 k1k2 . (2)
Equation ( 1) conserves the wave action N = n(k)dk and energy E = ω k n(k)dk. It has the zero flux (thermodynamic) solution n k = T /(k + µ) (with T and µ two constants), and the Kolmogorov-Zakharov isotropic constant flux solutions N (k) ∼ kn k ∼ k -2/3 and E(k) ∼ kω k n k ∼ k 0 which are the wave action and the energy isotropic spectra, respectively [6]. The correct dimensions can be recovered by noticing that the energy density follows in fact the relation [START_REF] Maggiore | Gravitational Waves[END_REF]

E ∼ c 4 32πG h 2 k 2 . ( 3
)
Complexity of equation ( 1) does not allow to derive analytically a simple model -a nonlinear diffusion equation -corresponding to the approximation of strongly local interactions. Below, we will assume that this type of interactions is dominant in GW turbulence and we shall derive two models by using phenomenological arguments.

III. FOURTH-ORDER DIFFUSION MODEL

Nonlinear diffusion models have proved to be a very useful tool in the analysis of both strong [10,16] and wave turbulence [12,13,19,[START_REF] Hasselmann | [END_REF]. Here, we shall derive such a model for GW turbulence. Since the leading nonlinear interaction of GW is four-wave interaction [6], the model is a fourth-order diffusion equation of the type

∂N (k) ∂t = ∂ 2 ∂k 2 D k N 4 (k) ∂ 2 (k 2 /N (k)) ∂k 2 , (4)
where N (k) is the one-dimensional (isotropic) wave action spectrum, k the wavenumber and D k a diffusion coefficient that eventually depends on the wave number k. This equation is constructed in such a way that it preserves the nonlinearity degree with respect to the spectrum (cubic in our case) and, its cascade and thermodynamic solutions. It is only in rare situations that one can derive a diffusion equation directly from the kinetic equation of wave turbulence (by taking the strongly local interactions limit). These exceptions concern nonlinear optics [11] and magnetohydrodynamics [13]. In our case, we will also use the phenomenology of wave turbulence [6]. A dimensional analysis of expression (4) gives the following information about

D k N (k) τ ∼ D k N 3 (k) k 2 , (5)
and thus

D k ∼ k 2 τ N 2 (k) . ( 6
)
We also have the scaling relation [6,[START_REF] Maggiore | Gravitational Waves[END_REF]]

E(k) ∼ h 2 k ∼ ωN (k) ∼ kN (k) , (7)
and therefore N (k) ∼ h 2 . This gives

D k ∼ k 2 τ h 4 ∼ k 2 (τ GW / 4 )h 4 ∼ k 2 τ GW ∼ k 3 , (8)
where the GW time is given by the relation τ GW ∼ 1/ω and where ∼ h µν /η µν ∼ h 1 is a small parameter. The introduction of expression ( 8) into (4) leads to the following fourth-order diffusion equation for isotropic three-dimensional GW turbulence

∂N (k) ∂t = A ∂ 2 ∂k 2 k 3 N 4 (k) ∂ 2 (k 2 /N (k)) ∂k 2 , ( 9
)
where A is a positive constant (a dimensionless number presumably of order one). Note that equation ( 9) is valid in the inertial range of weak GW turbulence. In the context of the primordial universe, the Planck length scale (∼ 1/k P ) can be seen as a small-scale limit for the diffusion model because the Einstein equations are not valid below, ie. when k > k P . On the other side, the largescale limit will be given by the scale k s at which GW turbulence becomes strong: it corresponds to space-time fluctuations h ∼ 1. Equation ( 9) conserves the wave action N (k)dk and energy ωN (k)dk. Indeed, if we define

∂N (k) ∂t = ∂ 2 K(k) ∂k 2 = - ∂Q(k) ∂k , (10)
with Q(k) being the wave action flux, then (for simplicity c = 1)

∂E(k) ∂t = k ∂N (k) ∂t = - ∂ ∂k K(k) -k ∂K(k) ∂k = - ∂P (k) ∂k , (11)
with P being the energy flux. We can check that the thermodynamic (zero flux) solution N (k) ∼ k 2 T /(k + µ) [6] is satisfied by equation ( 9). We can also find the constant (non-zero) flux solutions and find the cascade directions. Let us introduce N (k) = C N k α into equation ( 9); after simple calculations we obtain

Q(α) = -(1 -α)(2 -α)(3 + 3α)k 2+3α AC 3 N . (12)
Therefore, a constant wave action flux solution corresponds to α = -2/3. For this value we find Q 0 ≡ Q(-2/3) = -(40/9)AC 3 N < 0, which means that this solution corresponds to an inverse cascade (the wave action spectrum is positive definite thus C N > 0). Let us now substitute N (k) = C E k β into equation (11); after some calculations we obtain

P (β) = (1 -β)(2 -β)(-2 -3β)k 3+3β AC 3 E . ( 13
)
The constant energy flux solution corresponds to β = -1; in this case we have P 0 ≡ P (-1) = 6AC 3 E > 0, which means that we have a direct cascade of energy (because C E > 0). This analysis gives a prediction for the Kolmogorov constants C Q and C P , which depend, however, on A. We find

N (k) = (-Q 0 ) 1/3 C Q k -2/3 , (14)
with the Kolmogorov constant C Q = (9/(40A)) 1/3 and

E(k) = P 1/3 0 C P k 0 , (15)
with the Kolmogorov constant C P = (1/(6A)) 1/3 . Interestingly, the ratio of the Kolmogorov constants becomes independent of A:

C Q C P = 27 20 1/3 1.105 . ( 16
)
Note that a similar situation was also found with a differential model for two-dimensional hydrodynamic turbulence [18].

IV. SECOND-ORDER DIFFUSION MODEL

The nonlinear diffusion model of GW turbulence gets simplified if we do not include the thermodynamic solutions: it becomes a second-order diffusion equation which is easier to simulate numerically. We start with a diffusion equation of the type

∂N (k) ∂t = ∂ ∂k Dk N 2 (k) ∂(kN (k)) ∂k , ( 17
)
where Dk is a diffusion coefficient which may depend on k. As for the fourth-order diffusion equation, this equation is constructed in such a way that it preserves the nonlinearity degree with respect to the spectrum and its cascade solutions. A dimensional analysis of expression (17) gives

N (k) τ ∼ Dk N 3 (k) k , (18)
and thus (with

N (k) ∼ h 2 ) Dk ∼ k τ N 2 (k) ∼ k (τ GW / 4 )N 2 (k) ∼ k 2 . ( 19
)
Then, we obtain

∂N (k) ∂t = B ∂ ∂k k 2 N 2 (k) ∂(kN (k)) ∂k , (20)
where B is a positive constant (a dimensionless number presumably of order one). (See Section III for a discussion about the domain of validity of weak GW turbulence.)

We can check that the constant flux solutions discussed above are recovered by this equation, and that the wave action and the energy are conserved (the relation dP = kdQ can be used). By using similar notation as above, we can demonstrate that Q 0 ≡ Q(-2/3) = -(1/3)BC 3 N and P 0 ≡ P (-1) = (1/3)BC 3 E , which means that the directions of the cascades are recovered. For this model we obtain

N (k) = (-Q 0 ) 1/3 C Q k -2/3 , (21)
and

E(k) = P 1/3 0 C P k 0 , (22)
with equal Kolmogorov constants, C Q = C P = (3/B) 1/3 . Time-dependent solutions of equation ( 20) will be studied in section V. In particular we shall find a nonstationary solution with the wave action spectrum propagating towards small wavenumbers. This system is of finite capacity, ie. the integral

k I 0 N (k)dk , (23)
is finite when the solution ( 21) is used. This leads to an anomalous scaling with a non-trivial power-law. This non-stationary spectrum may be modelled as a selfsimilar solution of the second kind (see eg. [21,24]) taking the form

N (k) = 1 τ a N 0 k τ b , (24)
where τ = t * -t, and t * is a finite time at which the wave action spectrum reached the smallest available wavenumber. By introducing the above expression into (20) we find the condition

2a -b = 1 . (25)
A second condition can be found by assuming that N 0 (ξ) ∼ ξ m far behind the front. Then, the stationarity condition gives the following relation

a + mb = 0 . ( 26
)
Finally, the combination of both relations gives

m = - a b = - 1 2 - 1 2b . ( 27
)
The latter expression means that we have a direct relation between the power law index m of the spectrum and the law of the front propagation which follows k f ∼ τ b . For example, if we assume that the stationary solution -the Kolmogorov-Zakharov (KZ) spectrum -is established immediately during the front propagation [24], then m = -2/3 and b = 3 (and a = 2). In this case, the prediction for the front propagation is

k f ∼ (t * -t) 3 . ( 28
)
Any deviation from this prediction has to be considered as an anomalous scaling which is sometimes difficult to observe numerically because the deviation is often tiny [25,26].

V. NUMERICAL SIMULATION

In this section we shall investigate numerically the time evolution of the wave action spectrum described by the second-order diffusion equation ( 20) with B = 1. Linear hyper-viscous and hyper-hypoviscous terms are added in order to introduce sinks at small scale and large scale, respectively, to avoid the development of numerical instabilities. Then, the following equation is simulated

∂N (k) ∂t = ∂ ∂k k 2 N 2 (k) ∂(kN (k)) ∂k -νk 4 N (k) -η N (k) k 4 , (29)
with ν = 10 -85 and η = 10 20 . A logarithmic subdivision of the k-axis is used with k i = 2 i/10 and i an integer varying between 0 and 799. A Crank-Nicholson numerical scheme is implemented for the nonlinear term and an adaptive time-step is used. No forcing term is added. At the boundary cells a zero flux of wave action is imposed, then the (virtual) cells out of the numerical domain have a constant zero wave action. The code is publicly available from https://git.ias.u-psud.fr/ ebuchlin/nldiffus-gw.

Figure 1 shows the time evolution (every 6000 timesteps) of the wave action spectrum with an initial injection at k I ∼ 10 22 . As expected, an inverse cascade appears with a spectrum reaching the smallest scale available (k ∼ 1). An extended power law spectra (over more than 20 decades!) is formed with, however, a power law exponent slightly different from -2/3 as we can see at the bottom of Fig. 1 where compensated spectra are displayed.

The anomalous exponent can be determined very precisely by solving numerically an eigenvalue problem characterizing the existence of self-similar solutions of the second kind to the evolution equation. Such solutions were described in [21] for a broad class of second-order diffusion equations, exhibiting (forward) flux solutions with finite ultra-violet capacity. In the present case, a similar analysis carries through provided one maps the The black dashed line corresponds to (t * -t) 3.296 . inverse cascade behavior of the evolution equation ( 20) into a direct cascade in physical space. Explicitly, with the change of variables k → 1/ , N → Ñ = N/ 2 , the evolution equation (20) becomes

∂ Ñ ( ) ∂t = ∂ ∂ 4 Ñ 2 ( ) ∂( Ñ ( )) ∂ . (30)
This renormalized system has the equilibrium solution Ñ ∼ -1 and the flux solution Ñ ∼ -4/3 , which has finite ultra-violet capacity: the general framework of [21] then predicts the existence of a self-similar solution Ñ ( , t) = (t -t) -x/(2x-1) φ[ / ], where φ is a shape function describing the spectrum profile. The front evolves as (t) ∼ (t -t) 1/(2x-3) for a specific and uniquely determined exponent x = x , that characterizes the spectrum at small scales, Ñ ( , t) ∼ -x for 1. The anomalous exponent x is determined from a dichotomic search, a prescribed in [21]: We find x 1.3483, and this in turn implies the scaling exponent of the original system to be m = x -2 -0.6517. Quite remarkably, this is indeed the anomalous behavior observed in Fig. 2 (top) which shows the compensated spectra. Clearly the KZ solution is established only at t > t * , ie. after the largest (hyper-hypoviscous) scales are reached. This agrees with the time evolution of the wave action flux (bottom): a plateau corresponding to the KZ solution appears at the latest times only. Note that the simulation is stopped before the formation of the KZ spectrum over the entire range of scales.

We may also investigate the front propagation towards small wavenumbers and check if the power-law

k f ∼ (t * -t) 3.296 , (31)
corresponding to the anomalous scaling m -0.6517 (and also a -2.148 and b 3.296) is observed. For that, one needs to follow the front propagation k f (t) which will be defined by using the compensated spectra in Fig. 1 (bottom): the value 1.4 × 10 9 for the compensated spectra is chosen to define the front k f (t). The result (Fig. 3; inset) displays a sharp decrease of the wavenumber of the front at a time close to 10 -11 , which will be used to define t * . The value t * = 1.002168 × 10 -11 is taken because it corresponds to the most extended power law (dark blue curve). We see (Fig. 3) that the expected power-law (31) is well observed over six decades. This result illustrates the explosive character of the inverse cascade of wave action in GW turbulence. Note that the 24 decades in wavenumbers used for the simulation are necessary to detect without ambiguity the anomalous scaling. However, as explained above, it is expected that GW turbulence becomes strong at large scale, which prevents the formation of such an extended power-law range.

VI. CONCLUSION

In this paper we have proposed two nonlinear diffusion models for GW turbulence which reproduce the prop-erties previously derived from the kinetic equation [6] (power-law solutions, cascade directions). We have also derived some specific properties like the Kolmogorov constants. Additionally, we have performed a numerical simulation of the second-order diffusion model to illustrate the existence of an explosive inverse cascade of wave action, and have demonstrated the existence of an anomalous scaling which is typical to finite-capacity turbulence systems. This finding leads to a non-trivial power-law behind the propagating front in the inverse cascade spec-trum. This analysis of the anomalous spectrum in the case of an inverse cascade is the first made with a diffusion model. Note, however, that a similar analysis based on kinetic equations has already been done by [15,27,28] to study the formation of a Bose-Einstein condensate. The numerical simulation of the kinetic equation derived in [6] is much heavier to perform but not impossible. We can already predict the probable existence of an anomalous scaling for the wave action spectrum with, however, a possible difference in the power law index, as previously observed for other systems [14,21].

FIG. 1 .

 1 FIG.1. Time evolution of the wave action spectrum with (bottom) and without (top) a compensation by k 2/3 . The initial spectrum is localized at small scales (k ∼ 10 22 ). The comparison with the solution k -2/3 (black segments) reveals a slight difference.

28 FIG. 2 .

 282 FIG.2. Top: time evolution of the compensated (by k 0.6517 ) wave action spectrum; the Kolmogorov-Zakharov solution is established only after the hyper-hypoviscous scale is reached (the black segment corresponds to a spectrum in k -2/3 ). Bottom: time evolution of the wave action flux Q(k).

11 SlopeFIG. 3 .

 113 FIG.3. Temporal evolution of the spectral front k f for t ≤ t * in linear-logarithmic coordinates (inset) and as a function of t * -t in double logarithmic coordinates for three values of t * . The black dashed line corresponds to (t * -t)3.296 .