Jian Fu

Jean-Charles Maré

Yongling Fu

INCREMENTAL MODELING AND SIMULATION OF MECHANICAL POWER TRANSMISSION FOR MORE ELECTRIC AIRCRAFT FLIGHT CONTROL ELECTROMECHANICAL ACTUATION SYSTEM APPLICATION

Keywords: Bearing support translation, compliance contact, damping, elastic force [N]

In the field of more electric aircraft, electromechanical actuators (EMAs) are becoming more and more attractive because of their outstanding benefits of aircraft fuel reduction, maintenance costs saving, and system flexibility improvement. For aerospace electromechanical actuator applications, mechanical power transmission is critical for safety, in which reflected inertia to load, heat generated by energy losses and faults due to jamming, free-play and free-run are specific issues. According to the system-engineering process and simulation-aided design, this communication proposes an incremental approach for the virtual prototyping of EMA mechanical power transmission. Resorting to the Bond-graph formalism, the parasitic effects are progressively introduced and realism of models is increased step-by-step. Finally, the numerical implementations are presented and compared with basic, advanced and full models of mechanical power transmission in AMESim environment. Multi-level submodels are available and can be re-used for preliminary sizing, thermal balance verification and response to fault analysis.

NOMENCLATURE C b C m C j C s

Bearing translation support, motor electromagnetic, nut-screw inertial, rod output to surface torque [N/m] F 0

Preload force [N]

INTRODUCTION

Safer, cheaper and greener technologies are important initiatives for the next generation air transport in upcoming decades. In response to these needs, the aerospace industry is looking for an innovation (incremental or disruptive) in safetycritical actuation systems. In recent years, a significant interest is towards "more electric aircraft". The trend is to increase the usage of power-by-Wire (PbW) electrical actuators: electrohydrostatic actuator (EHA) and electro-mechanical actuator (EMA). These are envisioned to take the place of conventional hydraulic servo actuators (HSA). Compared to EHAs, EMAs totally remove the central and local hydraulic circuits, resulting in increased economic, competitive and environmental

IMECE2016-66436

advantages [START_REF] Botten | Flight control actuation technology for next-generation all-electric aircraft[END_REF]. EMAs have already served in a few secondary flight controls and for landing gear braking [START_REF] Todeschi | Airbus-EMAs for flight controls actuation system-perspectives[END_REF]. However, the maturity level of EMAs technology for primary flight control is far behind that of HSAs because of several potential issues: huge reflected inertia for dynamic performance, heat rejection for thermal balance and response to fault (jamming or freeplay) for failure. Therefore, EMAs still suffer from key issues that put them aside from being employed front line (normal or active mode) in safety-critical applications.

In this situation, for assessment of concepts during the preliminary design stages, as well as for performance virtual verification, system-level modelling and simulation has become well established. Despite the wide availability of commercial simulation software, it appears that there is a lack of methodology for the structuration of models. This particularly concerns the continuity between engineering activities, the knowledge capitalization, and the multi-physic nature of actuation. Therefore, it has become mandatory to compensate this lack by resorting to a model-based system engineering (MBSE) approach. This approach enables engineers with efficient means to access performance through virtual means and to increase maturity, performance and robustness. This communication applies to mechanical power transmission (MPT), which is a key part of EMAs. The designers and manufacturers must consider with great care the parasitic effects resulting from the imperfections of technology: friction, backlash, compliance, etc., and these considerations are especially true for safety-critical embedded power systems:

-consideration of energy losses for energy/power consumption and thermal balance;

-coupled and nonlinear effects, e.g. influence of temperature, influence of pre-loading on friction, influence of service on preload or backlash; -fault injections are considered for response to fault and fault-to-failure mechanism.

The following sections describe the power flows and multidomains disciplines in EMA system. By using systemengineering (SE) processes and simulation-driven design (SDD) [START_REF] Maré | Requirement-based system-level simulation of mechanical transmissions with special consideration of friction, backlash and preload[END_REF], MPT modeling is addressed to the current engineering activities and driven by requirements from engineering needs. Based on Bond-graph formalism, the replaceable, incremental and generic MPT models are developed by progressive levels of representativeness [START_REF] Fu | Incremental modelling and simulation of power drive electronics and motor for flight control electromechanical actuators application[END_REF]. Starting form a functional model that only considers the physical effects combined to perform the function, imperfections of technologies are progressively introduced as parasitic effects in the advanced and full MPT models. At last, these proposed models are implemented in causal simulation environment of AMESim and numerical simulations are compared and analyzed.

EMA SYSTEM DESCRIPTION

In this communication, a typical direct drive linear EMA is considered, in which the rotating nut of mechanical power transmission is directly integrated as part of the electric motor rotor. Due to elimination of the intermediate gear box, it is a compact design, and this "in-line" EMA is more attractive than geared EMA for aerospace application, thanks to the weight reduction and easy geometrical integration with-in the airframe. The EMA is developed for flight control application to drive a surface according to the pilot/autopilot commands, as described schematically in Fig. 1.

Motor

Sensor

Control surface

Position demand

Nut-screw

Motor

Power Drive Electronics

Velocity loop

Cycles rates of PWM The conventional architecture of an "in-line" EMA includes one electrical motor (EM) that converts the electric power into mechanical one, then the mechanical power transfers from rotation to translation domain through a nutscrew to move the control surface. Typically, the control of such an EMA is based on three main electronics devices. a) The flight control computers (FCC) ensure that the surface position corresponds to the flight law commands. Also for monitoring purpose, the treatment of faults monitoring from the actuator is completed by the FCC. b) The EMA actuator control electronics (ACE) receives the EMA position order from the FCC, sensor information from the EMA and processes an associated control loop. Usually, EMA control involves a rod position (outer loop), speed loop (middle loop) and a motor current/torque inner loop, for meeting the performance and stability requirements. ACE internally generates the pulse width modulation (PWM) command for the power drive electronics (PDE), as an image of the torque reference. c) The PDE is integrated near the EM in the EMA; modulated electric power forms the aircraft network to an EM.

Cross-linked Multi-domain Disciplines

Designing an EMA system requires multidisciplinary approaches not only for preliminary power sizing and estimation of the mass and geometrical envelope, but also for design consideration of thermal behavior. It is well known that the PDE produces conduction and switching losses in electrical domain, the EM produces copper loss and iron losses in electrical and magnetic domains and the MPT produces friction loss in mechanic domain. Unlike a HSA, where the fluid power also has a cooling effect, an EMA system is a thermal closed-circuit system. The heat generated by all mentioned power losses, is over dissipated towards the local environment. This may increase the EMA operating temperature and bring snowball effects within all the devices. By the above analysis based on Bond-graph theory [START_REF] Karnopp | System dynamics: modeling, simulation, and control of mechatronic systems[END_REF], the disciplines involved by EMA are steady in cross-linked multidomains, as shown in Fig. 2.

Electrical

EMA

Mechanic

Control

Thermal Magnetic

Functional Power Flow Analysis

In a direct drive linear EMA, there exists two functional types of motion on the same axis: rotation of the EM rotor, and translation of the MPT rod. An idea for virtual prototyping of EMAs is to make the same topology models as the sectional view of product geometrical arrangement. Figure 3 illustrates the proposed architecture and the corresponding functional power flow in an EMA model based on the Bond-graph formalism [START_REF] Karnopp | System dynamics: modeling, simulation, and control of mechatronic systems[END_REF]. The EM has been studied and presented in [START_REF] Fu | Incremental modelling and simulation of power drive electronics and motor for flight control electromechanical actuators application[END_REF]. The modelling of MPT is based on two degrees of freedom (2-DoF) approach [START_REF] Maré | 2-D lumped paramaters modelling of EMAs for advanced virtual prototyping[END_REF], but developed as an integrated model, in which the nut-screw, bearings (cylinder pairs), joints (hinges), shaft and end-stop are considered at the location of their power path.

Electric Motor (EM)

Housing with thermal model

Airframe

To ambiance

F V C  T Ṡ U m I m  Ṡ Rod To surface Fs Vs C  Ft Vt C 

Housing

Rod

Control

Mechanical Power Transmission (MPT)

F V C   Ṡ C m  m FIGURE 3. POWER FLOW IN A SIMPLIFIED TOPOLOGY STRUCTURE OF EMA.
In this concept of MPT, the thermal balance is considered, and the functions of rod anti-rotation and rotor axial thrust bearing are explicitly modeled. In the figure 3, for simplified thermal behavior study, EMA housing is assumed as the major and unique thermal body that receives the heat generated by power losses from EM and MPT and exchanges with the ambiance. The following achievements presented in this paper will focus on the MPT subsystem of the EMA.

Model Architectures versus Engineering Needs

In order to follow the system engineering (SE) process, the system-level EMA modeling and simulation should consider the physical effects at the right level of model complexity to meet the different engineering needs. From an engineering point of view, the best model of the technological device is usually not the most detailed one. In order to facilitate the numerical implementation and avoid hard nonlinearities and discontinuous effects, in some cases, the parasitic effects that weakly impact performance are neglected. In this manner, the above mentioned ideas are applied to the example of MPT. Table .1 shows that the MPT model architectures are driven by engineering needs and in support of the whole EMA system virtual prototyping.

Incremental Multi-level Modeling

The physical effects in MPT are mainly belonging to the mechanical domain, and some of them are coupled with nonlinear thermal behavior (e.g. temperature sensitivity). To address the requirement-based model development, an incremental modeling approach using the Bond-graph formalism is proposed, which step-by-step combines the parasitic effects (resistance, capacitance, inertance and etc.) and progressively increases the MPT model complexity (from functional to full, from linear to nonlinear, from continuous to discontinuous). All introduced multi-level models of MPT are scalable and replaceable (as far as possible) with common interfaces (power and signal).

 m C m V s F s FIGURE 4. SCHEMATIC OF A MPT MODEL (NUT-SCREW) [7].

BOND-GRAPH MODELING OF MPT

In direct drive linear EMAs, the MPT appears as a rotary to linear power transformer. This transformation is obtained, in this example, by using a key component: the nut-screw. It aims at converting the rotational velocity and torque of the EM into linear speed and force for driving the load, as shown in Fig. 4.

For lumped-parameter modeling of the MPT and according to the candidate location of the generic effects considered, the proposed decomposition involved five effects in series from motor to surface, as shown in Fig. 5: a global inertia, a perfect nut-screw, a friction loss, a compliance effect (that can represent backlash, preload, or pure compliance), and a global mass. In the proposed MPT architectures, friction and compliance effects are the most significant, and can be modeled at different levels of complexity (linear/nonlinear and continuous/discontinuous). Inertia and mass effects in the MPT can be optionally neglected or modeled depending on the order of the EM rotor and surface ones. The effects of the bearings, joints and end-stop will not be explicitly modeled but can be considered to make part of the MPT model. Consequently, the proposed generic model of the MPT involves four mechanical power ports (rotation and translation of nut and screw). An additional thermal power port enables the MPT to output the heat generated by its power losses. For advanced MPT models, the heat port provides the temperature that enables reproducing its influence on friction and compliance.

Perfect nut-screw Friction Compliance

Power losses for heat

v m  s v s F m F s C s  m C m Motor Surface Nut- support Screw- support MPT Inertia Mass FIGURE 5. PROPOSED ARCHITECTURE OF A MPT FOR BOND-GRAPH MODELING.

Functional Model

According to the above mentioned incremental approach, the MPT can be first expressed as a functional model with no parasitic effect. A TF element with causality in Bond-graph is shown in Fig. 6, in which the power flow is considered: The perfect nut-screw model achieves pure power transformation of ratio (2p/p), where p is the screw pitch.

TF p/2p  Perfect nut-screw
2 sm FC p p  (1)
2 sm p v  p  (2)
In the proposed causality, nut-screw receives the effort of torque (C m ) from EM and transfers the effort of force (F s ) as output to drive the surface.

Basic Model

This modelling level aims at increasing realism while being linear for control design activities. In this level of modeling, three major parasites are introduced: the screw inertia (J s ), the viscous friction loss (f e ), and the basic compliance effect including the contact elastic force (F e ) and damping force (F d ). Corresponding Bond-graph model is shown in Fig. 7, in which the proposed causalities keep the same interfaces of the power and signal ports as that of the previous functional model. Inertia Effect: J n is the screw inertia, I element in the Bond-graph, which generates an inertial torque:

n jn d C dt   (3)
Generally, the screw inertia, like the roller-screw type direct drive EMAs, is small enough when compared to the rotor/nut inertia. It can be neglected or integrated with respect to the rotor/nut as an overall inertia effect for MPT.

Friction Effect: In MPT, there are several solutions to model the friction loss: consideration of velocity, external load force and operating temperature sensitivity [START_REF] Maré | Friction modelling and simulation at system level: a practical view for the designer[END_REF]. From control engineer's point of view, the simplest one is to describe the friction loss that is proportional to the operating velocity, which is affected by a viscous coefficient (f e ), given by:

f v r F f v  (4)
Compliance Effect: Within EMAs, the MPT is not infinitely rigid. This makes it compliant due to the elastic deformation (x e ) of solids under mechanical stress, in particular at contact locations. A simple compliance model consists of a pure spring (k e ) and damping (d e ) effects. The contact elastic force (F e ) can be expressed as:

e e r F k x  (5)
Although the compliance damping is low, it must be considered in parallel with compliance in order to avoid unrealistic simulated oscillations. As the physical knowledge is very poor, the damping force (F d ) is usually considered as a linear function of the relative velocity (v e ), and of damping coefficient (d e ).

Advanced Model

This advanced level of MPT modeling is designed to offer high realism. In this manner, confidence is improved. For example, Energy balance studies are consistent from the energy balance point of view. Firstly, the nonlinear friction model (influence of velocity, load force and temperature) is upgraded to a modulated resistance element (MRS) in Bond-graph as it dissipates power. The nonlinear compliance model (pure spring, backlash and preloading effects), is considered as modulated capacitance element (MC). Then the mechanical power losses of friction (P f ) and compliance damping (P d ) can be introduced by the heat flows and connected to an additional thermal port. Consequently, it can easily enable the temperature to be used as a time variable input in the models of energy losses. Finally, the proposed MPT model architecture in this level involves a thermal port to be connected with EMA thermal model. The heat generation can be thought as an output to further EMAs' thermal model for whole system energy balancing. Thus the temperature sensitivities to the friction model and to the compliance model (dilatation) can be modelled using the temperature variable provided by the thermal port. Figure 8 describes the Bond-graph modeling of this advanced MPT model. Causalities are kept consistent with the two previous levels models.

Nonlinear Friction Effect:

The friction model can be built considering whether that function is velocity and load dependent, or load dependent and load independent [START_REF] Maré | Friction modelling and simulation at system level: Considerations to load and temperature effects[END_REF].

-Velocity and load dependent: for nut-screws, a five parameter model, Eqn. [START_REF] Cochoy | Integration and control of electromechanical and electrohydraulic actuators in a hybrid primary flight control architecture[END_REF], has been identified by Karam [START_REF] Karam | Modelling and simulation of mechanical transmission in roller-screw electromechanical actuators[END_REF]. The friction force introduces a constant Coulomb friction (first part), a Stribeck effect of rapid change of friction at low velocity (second part), and a load and power quadrant dependent Coulomb effect (last part):

( sgn( ) sgn( )

r st vv f cl st L L r r F F F e F a b F v v        (7)
-Load dependent and load independent: this model is consistent with nut-screw suppliers datasheets, which provide efficiency (the load dependent friction), as well as the no-load friction under opposite load and the no-drive friction under aiding load (the load independent friction). The velocity effect is added in a second step by introducing its influence of these parameters. The details of this modeling approach has been presented by author in [START_REF] Maré | Friction modelling and simulation at system level: Considerations to load and temperature effects[END_REF].

For any friction model in MPT, the power loss (P f ) is calculated as

f f r P F v  (8)
Nonlinear Compliance Effect: From a mechanical engineer's point of view, obtaining a realistic model of compliance is of particular importance because compliance can significantly influence the low-magnitude dynamic performance [START_REF] Fu | Preliminary design rules for electromechanical actuation systems -effects of saturation and compliances[END_REF] and the service life of the EMAs. Backlash or preload are two key effects in a realistic compliance model. In the absence of preload and backlash, the axial mechanical stiffness will suffer from such a reduced stiffness around the null transmitted force where not all the contacts are fully loaded, this generates the so-called "lost motion". Figure.9 shows the elastic characteristics of proposed compliance model (either pure spring, backlash or preload) and introduced by a single parameter x 0 : -Pure spring effect: used when x 0 = 0, the contact force F c is purely proportional to the relative displacement x r

            (11)
For numerical stability and rapidity, it is proposed to implement the compliance model by combining Eqns.( 9) to [START_REF] Fu | Preliminary design rules for electromechanical actuation systems -effects of saturation and compliances[END_REF] in such a way so as to avoid switches or "if" functions:

0 0 0 (1 ( ) ( ) 2 ()1 ( ) 2
                 (12)
Most of the common contact models implemented in simulation SW fail to avoid discontinuity of contact force when the contact is reached or lost. This is due to the damping effect that is not implemented with care (damping force is null on reaching the contact and contact force is never attractive). Thus the damping force F d can be modeled in an advanced manner as a force which is acting in parallel with the elastic force under the following constraints:

F dv x x           (13)
The power loss due to compliance damping (P d ), although negligible in general, can be calculated to make the MPT model exactly balanced with respect to energy:

d de P Fv  ( 14
)
Temperature Sensitivity: In MPT, temperature influences friction due to viscosity effects and dimensions due to dilation.

-Influence on friction: as it is well known, most EMAs' MPT are grease lubricated. High temperatures lowers the viscosity of lubricant, resulting in increment of the friction loss. However, the effect of temperature on friction is poorly documented in the nut-screw suppliers' datasheets. At systemlevel modeling, several models have been proposed by author in [START_REF] Maré | Friction modelling and simulation at system level: Considerations to load and temperature effects[END_REF]. As a consequence, a linear or nonlinear parameter/function ( f  a consequence, a linear or nonlinear ) can be introduced to modulate the friction force: ( ( sgn( ) sgn( ))

r st v f f cl st L Lr r F F Fe F a b Fv v         f f cl st F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st f f cl st f f cl st f f cl st F F F e F a b F v v f f cl st f f cl st F F F e F a b F v v f f cl st F F F F a b F v v      F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st      f f cl st F F e F a b F v v f f cl st F F F e a b F v v F F F F a b F v v    F F F e F a b F v v F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st f f cl st F F F e F a b v f f cl st      f f cl st F F e F a b F v f f cl st f f cl st F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st f f cl st F F F e F a b F v v f f cl st F F F e F a b F v v      F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st      f f cl st F F e F a b F v v f f cl st F F F e F a b F v v F F F e F a b F v v f f cl st f f cl st f f cl st F F F e F a b F v v f f cl st f f cl st F F F e F a b F v v f f cl st f f cl st f f cl st f f cl st F F e F a b F v v f f cl st f f st F F F e F a b F v v f f cl st F F F e F b F v v      F F F e F a b F v v F F F e F a b F v v      F F F e F a b F v v f f cl st F F F e F a b F v v f f cl st     f f cl st F F e F a b F v v f f cl st f f cl st F F F e F a b F v v f f cl st      f f cl st F F F e F a b F v v f f cl st F F F e F a b F v v F F F e F a b F v v     F F F e F a b F v v F F e F a b F v v f f cl st F F F e F a b F v v f f cl st f f cl st F F F e F a b v v f f cl st      f f st F F F e F a b v v f f cl st f f cl st F F e F a b F v v f f cl st ( 15
)
-Influence on dimension: in MPT, especially for roller type, the nut-screw has a very small lead. The thermal expansion may cause the thermal dilatation. This leads to variation of dimensions and may impact backlash, preload and friction. Once again, if a model of dilation is developed, this effect is easily introduced by using the temperature variable of the MPT thermal port. For example, the temperature effect on the value of the backlash parameter x 0 can be a linear or nonlinear parameter/function ( x  of the backlash parameter ), and to describe the relative dilatation effect, following formula can be adopted: 0 0

(1 )

x x     0 0 (1 ) 0 0
(1 )

0 0 0 0 0 0 0 0 (1 ) 0 0 0 0 (1 ) 0 0 0 0 0 0 0 0 0 0 0 0 (1 ) 0 0 0 0 (1 ) 0 0     0 0     0 0 (1 )     (1 ) 0 0
(1 )

0 0    0 0
(1 )

0 0 0 0 (1 ) 0 0 0 0 (1 ) 0 0     0 0
(1 )

0 0 0 0 (1 ) 0 (1 ) (1 )     (1 ) (1 ) 0 0
(1 )

0 0 0 0 (1 ) 0 0     0 0
(1 )

0 0 0 0 ) 0 0 (16)

Full Model

This full model available involves mechanical faults (e.g. jamming and free run) and parasitic motion (relative displacement because of bearings and joints), shown in Fig. 10. Mechanical Faults: It is definitely an added advantage when the model can reproduce the major faults that may occur in MPT. Jamming in MPT may jam the entire control surface. It can be modeled as a signal input to modulate the friction model parameters (MRS). If wear of the mechanical components will cause control surface free-play, this failure can be also modeled in a manner as signal modulates the model of compliance (MC): preload is reduced and backlash is increased by modifying parameter x 0 .

-Jamming to failure: it can be modeled as brake function by adding a stiction force to the friction model:

* f jam f F F F  (17)
-Free-play or wear to failure: the compliance parameter x 0 can be varied in the presence of increased backlash or decreased preload, which depends on the fault signal input:

* 0 0 w x x x   (18)
Parasitic Motion: Another highlight in this level of MPT's architecture enables better modelling of mechanical interfaces: anti-translation function of motor/nut and anti-rotation function of rod/screw. In this way, it is possible to consider imperfect bearings and joints (e.g. compliance and friction) and to access the reaction torque (C b ) and force (F b ) to EMA housing or airframe, and the relative rotational velocity ( nr ) for nut/rotor and relative velocity (v sr ) for screw to surface can be modeled.

-Support rotation: the relative rotational velocity between the rotor/nut and support bearings can be modeled as:

nr n b      (19)
-Support translation: the relative translational velocity between the screw and support bearings given:

sr s b v v v  (20)
All these advances proposed finally result in the two degrees of motion (2-DoF) MPT model, which can connect the bearings support interfaces.

IMPLEMENTATION OF MPT MODELS

For illustrating the interest of the proposed MPT model, the flight surface is subjected to the aerodynamic load. The control of the EMA is based on cascade multi-loops control architecture, shown in Fig 11 . The PDE and EM are assumed to be perfect, the parasitic effects are neglected. Therefore, the PDE and EM are considered as a modulated source of effort: a second order function for PDE considering current/torque loop dynamics and a source effort of torque (S e ) for the EM based on Bond-graph modeling [START_REF] Fu | Incremental modelling and simulation of power drive electronics and motor for flight control electromechanical actuators application[END_REF].

The proposed MPT multi-level models are independent, and replaceable (same input power port from the EM, and output power port to the surface). Although causal choices have been made, each architecture of MPT is adapted in causal (AMESim) and non-causal (Dymola) simulation environment.

EMA Controller

Basic MPT Model

The basic model is always simple, low order and linear. Its use is relevant for basic simulation and for preliminary control synthesis. According to Bond-graph model in Fig. 7, the implementation of this basic MPT model in AMESim is shown in Fig. 12. In this basic MPT model, all the components come from standard libraries. The motor rotor and nut-screw inertia are merged into a single inertia on which applies a pure viscous friction. The motor power flows to a perfect nut-screw and is transmitted to the load through a compliance effect that is considered using linear spring-damper model. At MPT interface with load, a rod mass is also considered.

Advanced MPT Model

According to the Bond-graph model in Fig. 8, firstly, the former basic model of linear friction and compliance effects are replaced by more realistic models, Fig. 13. In advanced MPT model, the proposed friction model of Eqn ( 7) is implemented by making the friction force dependent on transmitted force, temperature and sliding velocity (force and velocity are captured through sensors models). This dependence can be either defined by parametric functions or using loop-up tables. The compliance model of Eqn. ( 12) is implemented by moving to the "signal" world, separating explicitly the flow and effort variables though a standard AMESim model. Secondly, a standard AMESim thermal power bond is introduced. It collects the power losses (friction model and structural damping of compliance model) and inputs temperature to the friction (impact on friction factor) and compliance (impact on backlash) models.

Full MPT Model

Based on the Bond-graph model in Fig. 9, the implementation of full MPT model in AMESim is shown in Fig. 14. The jamming and free-run faults are optionally introduced by external signal that affect the friction and compliance parameters, Eqns. ( 17) and (18), respectively. In addition, the relative rotational and translation positions between two mechanical components are available. Two mechanical ports associated with nut-housing for antitranslation and with screw-housing for anti-rotation, can be explicitly introduced for connection with the bearings and joints for 2-DoF modeling in future study.

Nut housing

Perfect nut-screw

Jamming

Free-run

To thermal model Source: i-product data [START_REF] Exlar | [END_REF],ii-scaled from former experiments [START_REF] Wang | A force equalization controller for active/active redundant actuation system involving servo-hydraulic and electro-mechanical technologies[END_REF], iii-literature [START_REF] Karam | Modelling and simulation of mechanical transmission in roller-screw electromechanical actuators[END_REF] Simulation parameters are shown in Tab.2 and Tab.3. The flight control surface is simply modelled as an equivalent translating mass (M s ) to which the air load is applied. In EMA controller, the speed limitation ( lim ) and torque limitation (C lim ) are introduced. G p (s) is the position controller and G v (s) is the velocity controller. In addition, the anchorage of the EMA housing to the wing and the EMA rod to load connection is all assumed to be rigid.

Interest for Control Design

An aileron position step demand of 10 mm is applied at time t=0.1 s, then followed by a step aerodynamic disturbance force of 10 kN at time t = 1 s. Figure 15 compares the load position simulated by the different MPT models (functional, basic and advance). Both simulated responses are stable. It can be seen that the realistic friction effect seriously increases EMA system damping and affects dynamic performances. Figure 16 displays that the motor torque/current is saturated for a longer time when a more realistic friction model is considered. Significant differences in response also appear at very low velocity.

Interest of Power Consumption Analysis

The proposed full model of MPT can be used for comparative analysis of the power losses and energy consumption for whole EMA system. Figure 17 (A) shows the mission (X c and F ex ) and the load position response (X s ): A trapezoidal position profile is demanded that of a maximum speed of 125mm/s. The external force is increased from 0 N at 0.1 s to 15 kN at t = 1 s. Focusing on MPT frictional loss, the losses of PDE and EM are adopted from advance model in Ref. [START_REF] Fu | Incremental modelling and simulation of power drive electronics and motor for flight control electromechanical actuators application[END_REF]. It can be seen from Figure 17 (B) that when the final surface position is reached (no speed but high load), friction loss is null; high speed and high adding load has the highest power losses. The MPT friction loss represents the highest source (70%) of total energy losses, which highlights the importance to developan advanced friction model of MPT.

Interest for Wear/Ageing and Response to Faults

The analysis illustrates the interest of the proposed models for wear/ageing (increased backlash), and free-run or jamming faults. Figure 18 (A) compares the simulated responses of the MPT full model, and introduces the functions: null preload/backlash, (x 0 = 0), a backlash (x 0 = 0.1 mm), then a preload of 3 kN (x 0 = -0.06 mm). It is obviously observed that

F f FF

 f Initial normal, faults injections and temperature sensitivity friction force [N] F L Load force [N] F m F s Motor shaft output, rod output to surface force [N] jam Jamming stiction force [N] f v Viscous friction coefficient [N/(m/s)] I m I s Motor windings, DC supplied current [A] P d P f Damping, friction loss [W] S Heat power [J/sK] U m U s Motor wingdings, DC supplied voltage [V] v b v e v r v s v sr Bearing support, elastic, relative, rod output to surface, relative rod/support translational velocity [m/s] Initial normal, faults injections, temperature sensitivity backlash/preload parameter [m] x c Position command [m] x w Wear parameter of nut-screw [m] x r Relative elastic deformation [m]  f f  Initial, temperature sensitivity friction factor[-] x   Temperature dependency backlash/preload parameter [-]  d  i Efficiency direct, indirect [-]  Temperature [°C]  b  m  n  sr Support rotational, motor rotor, relative nut/support, relative rod/support angular velocity [rad/s]

FIGURE 6 .

 6 FIGURE 6. FUNCTIONAL MODEL OF MPT IN BOND-GRAPH.

FIGURE 7 .

 7 FIGURE 7. BASIC MODEL OF MPT IN BOND-GRAPH.

FIGURE 8 .

 8 FIGURE 8. ADVANCED MODEL OF MPT IN BOND-GRAPH.

FIGURE 9 .

 9 FIGURE 9. ELASTIC CHARACTERISTIC OF PROPOSED COMPLIANCE MODEL.

--

  Backlash effect: used when x 0 > 0, model that displays a total dead-zone of 2x 0 and the contact elastic force (F e ) Preload effect: used when x 0 < 0, where the preload force is 0

FIGURE 10 .

 10 FIGURE 10. FULL MODEL OF MPT IN BOND-GRAPH.

FIGURE 11 .

 11 FIGURE 11. SCHEMATIC DIAGRAM OF EMA SYSTEM FOR MPT MODELS IMPLEMENTATION

(FIGURE 12 .

 12 FIGURE 12. IMPLEMENTATION OF BASIC MPT MODEL

FIGURE 13 .

 13 FIGURE 13. IMPLEMENTATION OF ADVANCED MPT MODEL

FIGURE 14 .

 14 FIGURE 14. IMPLEMENTATION OF FULL MPT MODELNUMERICAL SIMULATION AND RESULTS

ForceFIGURE 16 .

 16 FIGURE 15. DYNAMIC PERFORMANCES COMPARISON.

  PDE losses (conduction and switching) EM copper losses EM iron losses (eddy current and hystesis) MPT losses (friction and damping) POWER CONSUMPTION AND LOSSES ANALYSIS FIGURE 17. EMA TOTAL POWER LOSSES (ADVANCED MODEL).

Proceedings of the ASME 2016 International Mechanical Engineering Congress and Exposition IMECE2016 November 11-17, 2016, Phoenix, Arizona, USA

TABLE 1 .

 1 MODEL ARCHITECTURES VS ENGINEERING NEEDS.

					Engineering needs
	Model architectures	No parasite	Power sizing	Thermal balance	Natural dynamic	Stability accuracy	Consumed energy	Failure response	Load propagation	Reliability
	Functional model
	perfect transformer (TF)	Y	Y	Y	Y	Y	Y	Y	Y	Y
	Basic model
	inertia/mass (I)	N/A	Y	N/A	P	Y	Y	N/A	P	P
	simple friction (RS)	N/A	P	N/A	P	Y	Y	N/A	P	P
	simple compliance (C)	N/A	P	N/A	P	Y	P	N/A	P	P
	Advanced model
	modulate friction (MRS)	N/A	Y	Y	Y	Y	Y	Y	Y	P
	modulate compliance (MC)	N/A	Y	Y	Y	Y	Y	Y	Y	P
	Full model
	failures injection	N/A	P	N/A	Y	Y	N/A	Y	Y	Y
	support interfaces	N/A	P	N/A	N/A	P	N/A	N/A	Y	P

Note: Y means yes; P means possible, but depends on relative level; N/A means not applicable

TABLE 2 .

 2 EMA CONTROLLER, PDE, LOAD PARAMETERS

	Symbol	Characteristic	Value
	EMA controller
	K p	Position loop proportional gain (radmm -1 /s)	45
	K v	Velocity loop proportional gain (Nms/rad)	0.47
	 lim	Velocity limitation (rad/s)	314
	C lim	Torque limitation (Nm)	10
	PDE dynamics parameters
	 i	Current/torque loop damping factor (-)	0.7
	f i	Current/torque loop natural frequency (Hz)	600
	Load parameters
	M s	Equivalent flight surface mass (kg)	600

TABLE 3 .

 3 MPT PARAMETERS FOR FULL MODEL

	Symbol Characteristic	Value
	Inertia effect i
	J n	Nut-screw inertia integrated with rotor (kgm 2 )	0.0017
	M t	Rod mass (kg)	1
	Perfect nut-screw i
	p	Lead of screw (mm)	2.54
	Friction effect ii
	F cl	Coulomb friction force (N)	7590
	F st	Stribeck friction force (N)	-4702
	v st	Reference speed for Stribeck friction (m/s)	0.035
	a	Mean coefficient of external force (-)	0.218
	b	Quadrant coefficient (-)	-0.13
	Compliance effect iii
	k e	Elastic stiffness of nut-screw (N/m)	310 8
	d e	Damping factor of nut-screw (Ns/m)	110 4

the backlash makes surface chattering and impacts the position accuracy. The preload removes this effect (but increases friction) at the surface and affects the surface rapidity.

CONCLUSIONS

An incremental modelling and simulation of mechanical power transmission has been proposed for supporting virtual prototyping of EMAs with a system-level view of MBSE design. Combining perfect nut-screw, parasitic effects (friction and compliance) and fault failures; the functional, basic, advanced and full models have been developed from the needs generated by engineering tasks and intends to use objectoriented elements and interfaces from standard model libraries. Each lumped-parameter model is energy balanced and mechanically replaceable (interface). It enables faults (jamming or free-play) to be simulated. The models are ready to be extended for next step thermal simulation and to be integrated in full incremental realistic more electric aircraft virtual prototyping.


Karin Müller

Céline Duwig

Anne-Julie Tinet
email: anne-julie.tinet@univ-lorraine.fr

Alfonso Gastelum

Lorenzo Spadini

Marie Christine Morel

Pascal Charrier

Orchard management and preferential flow in Andosolscomparing two kiwifruit orchards in New Zealand

Keywords: carbon, filtering function, kiwifruit orchard, soil architecture, X-ray computed tomography

Sustainable horticulture depends on the integrity of soil functions, which directly depend on soil architecture affecting aggregation, root growth, as well as liquid and gas permeability. We hypothesised that changes in soil architecture resulting from feedback mechanisms between management, soil organic carbon contents (SOC), biota and vegetation can be captured with X-ray computed tomography (CT), and that these affect the soil filtering function, which thus, can be manipulated through orchard management. We compared the transport of copper, a widely used fungicide, through intact soil cores from vine rows of kiwifruit orchards under organic and integrated management. We first derived 3D-macropore characteristics from CT-images, followed by leaching a pulse of copper and a tracer through the same cores. The organic orchard soil had a significantly higher SOC content than the integrated orchard soil, and this was positively correlated with total porosity. Macropores (>92 µm) were larger with a higher connectivity, but significantly fewer in the organic than the integrated orchard soil. This resulted in a lower macroporosity and a better copper filtering capacity of the organic than the integrated orchard soil. Copper distribution was reasonably predicted when combining SOC contents, pH and macropore characteristics. Significant relationships between soil parameters and indicators of the strength of preferential flow verified that CT-derived macropore characteristics can be used to predict functional solute transport parameters. The relevance of our results and relationships observed between macropore characteristics, functional indicators of preferential flow and the fate of copper needs verification with samples representing more soils and sites.

Introduction

Kiwifruit is the fruit of a perennial woody vine. New Zealand (NZ) is among the leading kiwifruitproducing countries worldwide. Kiwifruit production totalled ~1 412 351 MT in 2012 with Italy and NZ together producing ~54% of this total (http://faostat.fao.org/site/339/default.aspx). In NZ in 2016, 9336 ha of Actinidia deliciosa 'Hayward' were under integrated management, which can be understood as a production system aiming at 'environmental and economic sustainability including the efficient use of natural resources' (http://www.zespri.com/), and 537 ha were grown according to certified organic standards (BIO-GRO; NZ's organic protocol organisation;

https://www.biogro.co.nz/).

Copper (Cu) has been applied to various crops for controlling different plant diseases since 1885, particularly in vineyards and orchards [START_REF] Alloway | Heavy metals in soils -trace metals and metalloids in soils and their bioavailability[END_REF]. For example, Bordeaux mixture, a combination of copper sulfate and lime, has traditionally been used for the control of fungal diseases in grapes. It is commonly accepted that long-term foliar applications of copper-based fungicides and bactericides can enhance Cu concentrations in orchard soils [START_REF] Nagajyoti | Heavy metals, occurrence and toxicity for plants: a review[END_REF]. Many studies have investigated Cu accumulation in orchard soils worldwide [START_REF] Besnard | Influence of organic amendments on copper distribution among particlesize and density fractions in Champagne vineyard soils[END_REF][START_REF] Pietrzak | Copper accumulation, distribution and fractionation in vineyard soils of Victoria, Australia[END_REF][START_REF] Vavoulidou | Copper content in agricultural soils related to cropping systems in different regions of Greece[END_REF][START_REF] Sun | Investigation of the copper contents in vineyard soil, grape must and wine and the relationship among them in the Huaizhuo Basin Region, China: a preliminary study[END_REF]. Background concentrations of Cu in soil are dependent on the soil parent material and are generally below 30 mg kg-1 soil [START_REF] Alloway | Heavy metals in soils -trace metals and metalloids in soils and their bioavailability[END_REF]. [START_REF] Besnard | Influence of organic amendments on copper distribution among particlesize and density fractions in Champagne vineyard soils[END_REF] reported 100-1500 mg Cu kg-1 in Champagne vineyard soils compared with 5-30 mg Cu kg-1 in adjacent agricultural soils. Horticultural soils with levels exceeding 3000 mg Cu kg-1 have been reported in a >100-year-old Brazilian vineyard [START_REF] Mirlean | Metal contamination of vineyard soils in wet subtropics (southern Brazil)[END_REF].

The motivation for this research was the significant increase in the use of antimicrobial Cu sprays on the vines in NZ kiwifruit orchards as a response to the first outbreak of the bacterial disease Pseudomonas syringae pv. actinidiae (Psa) in 2010 that has been threatening the NZ kiwifruit industry [START_REF] Vanneste | Recent advances in the characterisation and control of Pseudomonas syringae pv. actinidiae, the causal agent of bacterial canker on kiwifruit[END_REF]. 'BIO-GRO' orchardists in NZ are currently permitted to use 6 kg ha-1 year-1 of Cu, which is 3 kg more than the maximum allowed before the incursion of Psa (https://www.biogro.co.nz/). New Zealand's conventional growers producing for the international market are permitted to use 8 kg ha-1 year-1 of Cu. Several recent monitoring studies on NZ kiwifruit orchard soils have reported elevated Cu levels [START_REF] Dean | Effects of copper sprays on microbial communities in kiwifruit orchard soils[END_REF]. Increased Cu concentrations in soils enhance the risk of Cu leaching, which has been shown to be closely related to soil structure [START_REF] Paradelo | Effects of past copper contamination and soil structure on copper leaching from soil[END_REF]. For example, [START_REF] Paradelo | Effects of past copper contamination and soil structure on copper leaching from soil[END_REF] conducted displacement experiments with intact soil cores collected from a long-term contaminated pasture site with a Cu gradient of 20-3800 mg Cu kg-1 soil. They reported a strong correlation between Cu leaching and soil macroporosity derived from X-ray computed tomography (CT) scans with a resolution of 200 µm.

It is well known that carbon (C) contents affect soil structure and vice versa. The C contents in agricultural soils can be manipulated through soil management practices [START_REF] Deurer | The impact of soil carbon management on soil macropore structure: a comparison of two apple orchard systems in New Zealand[END_REF]). The link between soil biotic activity, soil organic matter decomposition and stabilisation, and soil aggregate dynamics has been recognised and intensively studied since the beginning of the 1900s [START_REF] Six | A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics[END_REF]). Soil porous structure and organic matter are not independent. Pores are created by abiotic and biotic factors (root growth and burrowing etc.) and different forms of organic matter stabilise pores of different sizes. In turn, pore characteristics influence organic matter dynamics [START_REF] Kay | Conservation tillage and depth stratification of porosity and soil organic matter[END_REF].

Applying 3D X-ray CT, it was found that soil organic carbon (SOC) and organic fertilisations influenced soil structure both in terms of pore morphology and pore network topology [START_REF] Ferro | Soil macro-and microstructure as affected by different tillage systems and their effects on maize root growth[END_REF]. [START_REF] Naveed | Impact of longterm fertilization practice on soil structure evolution[END_REF] evaluated the soil water characteristics and gas transport and analysed soil structure via X-ray CT. Their results suggested that with increasing fertiliser application rates, the pore size distributions widened, and pore connectivity was significantly improved. [START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF] found that preferential transport in macropores was reduced in soils with larger contents of SOC. These soils were also characterised by larger volumes of pores <600 μm in size.

We hypothesised that changes in pore structure resulting from the feedback mechanisms between orchard management, soil biota and vegetation can be captured with 3D X-ray CT, and that such changes can explain differences in the filtering function of soils. Accordingly, our objectives were (i)

to assess the effect of orchard management on topsoil C contents and architecture, and (ii) to analyse the relationship between the soil macropore network, water flow and Cu accumulation in Andosols, a typical soil for kiwifruit production in NZ.

Materials and methods

Sites and soils

The two selected sites are kiwifruit orchards developed on Andosols (IUSS Working Group WRB, 2006) located in the Bay of Plenty, NZ. The orchards are located in the same region (Orchard 1,176.05841;Orchard 2,176.40095) sharing the same climatic conditions.

Long-term average annual temperature in the region is 14.7°C and annual rainfall totals on average 1460 mm. The annual potential evapotranspiration rate in Te Puke is 1014 mm (Te Puke, NIWA climate station). The kiwifruit vines in the two orchards are of the same cultivar Actinidia deliciosa 'Hayward' and of similar age (29 and 25-30 years). One of the two orchards (O) has been under certified organic management (BIO-GRO) since its establishment. The second orchard (I) has been managed following integrated management principles, e.g. the 'KiwiGreen' program, an integrated pest management program, that ensures that the environmentally safest available pesticide is applied upon monitoring of occurrence or level of a particular pest or disease as well as weather conditions for its potential application. The soils developed from the same parent material. A preliminary analysis of generic soil properties of the two soil profiles to a depth of 1 m is summarised in Table 1, which confirmed that both soils had a silt loam texture. It also confirmed our expectation of significantly higher SOC contents and lower bulk densities in O than I.

The main differences between the two orchard management practices relevant to this project relate to those affecting SOC contents and thus in turn soil structure: in O ~6 t year-1 of an organic soil amendment consisting of 50% compost and 50% chicken manure have been applied since the establishment of the orchard, while in I this practice has started only in recent years (3-5 years) and at a lower rate. BIO-GRO regulations inhibit the use of herbicides. Accordingly, O's orchard floor was under full grass cover. In contrast, the vine rows in I were kept weed-free by annual applications of a glyphosate-containing herbicide at the maximum recommended application rate (660 g a.i. ha-1).

Intact soil cores from the topsoil of the vine rows were collected on 5th March 2012. A cylindrical PVC tube with an inner diameter of 104 mm and a length of 111 mm was pushed a few centimetres into the ground by hand until the penetration resistance was too high. The remainder of the core was carefully excavated by removing soil around the PVC tube, which allowed pushing it into the soil by hand. We also took bulk soil near each core to determine its gravimetric water content. The soil cores were weighed at field water content. We calculated each core's dry soil weight, bulk density and total porosity taking into account the PVC tube's tare weight and the soil water content at the time of sampling and assuming a particle density of 2.65 g cm-3. On 10th April another set of intact soil cores were collected to determine the initial total Cu concentration in the topsoils. All soil cores were stored at 4°C in the dark before the experiments.

X-ray CT scanning and image analysis

The scanning of the nearly saturated soil cores with 3D X-ray CT was conducted at the 3S-R laboratory in Grenoble, France. The scanning was performed with a multi-energy X-ray source at 120 kV and 500 µA. In each scan 1200 angular projection images were collected and each radiograph was averaged over 10 frames. Ring artefacts were minimised during data acquisition and a 0.5 mm Cu filter was used to reduce beam hardening. The software package X-act (RX-solutions, Hurst, Texas, United States) was used to obtain the 3D maps of attenuation coefficients from the 2D angular projections. The 3D images of attenuation coefficients with the isotropic voxel size of 92 µm were then translated into a continuous stack of 2D 16-bit TIFF images using the software DigiCT 2.42 (Digisens 3D, Le Bourget-du-Lac, France). The algorithm FDK was used to reconstruct the slices.

Voxel size was 92.94 × 92.94 × 92.94 µm. The images were binarised and the individual pores labelled using the module fuzzycmeans clustering method of the program SMAS (Gastelum-Strozzi 2018). In brief, the following four values were used to label void voxels: three texture values, namely the intensity value, the average intensity value, the homogeneity with respect to the neighbouring voxels, and the intensity value scaled with respect to the beam hardening. A pore was defined as a set of voxels with a six-connectivity, meaning that the neighbouring voxels in the three main directions within an image were considered. To derive morphological descriptors important for flow processes, like tortuosity and curvature, a medial axis representation of pores was defined as thin segment of connected voxels equidistant to its boundaries (skeletonisation, [START_REF] Ballard | Representation of three-dimensional structures[END_REF].

Only pores that could be defined by a skeleton were included in the analysis. A total of 14 morphological parameters were obtained for each pore of each core using the software SMAS (www.biocomlab.com/apps/SMAS/). In the following some of the parameters are explained. The circularity C of a pore was determined as:

2 4 A C P   (1)
where A is the area of a pore in the plane (m 2 ) and P is the perimeter (m). The tortuosity T is defined as the ratio between the geodesic length of the skeleton and the Euclidian distances between the highest skeleton vertex with respect to the z-axis and the lowest one. The Euclidian distance (m) is defined as:

  2 μμ E s s s s      (2)
where 𝑠 is the pore skeleton, and 𝜇 is the highest skeleton vertex and 𝜈 the lowest. The geodesic length (m) can be calculated as:

1 1 s iN ii i L s s      (3)
where 𝑁 𝑠 is the total number of vertices in the skeleton 𝑠. The tortuosity is then obtained as:

L T E  (4)
The connectivity was determined using the plugin MorphoLibJ (http://imagej.net/MorphoLibJ; [START_REF] Legland | MorphoLibJ: integrated library and plugins for mathematical morphology with ImageJ[END_REF] of the software ImageJ. The Euler number, which may be associated with the number of poresnumber of loops + number of solid inclusions, was divided by the total soil volume (V) to give a connectivity density CE . To obtain parameters for entire soil cores or segments of the core, the pore parameter values were averaged for each core or the respective segments by the number of slices. From here onward the pores detectable in the X-ray CT scans are referred to as macropores and they constitute the soil macroporosity. We defined percolating pores as those connected from the top of a core to any of the boundaries of the core.

Transport experiments

After completion of the scanning, we conducted transport experiments with the same soil cores close to saturated conditions. We used KBr as the water tracer and CuN2O6 as the solute. The integrity of the soil structure was ensured by placing a grid at the bottom of the cores. A rainfall simulator was placed on top of the cores. It consisted of nine hypodermic needles fed from a reservoir.

Two peristaltic pumps (Gilson Minipuls 3, Gilson, US) were used to create steady-state flow through the soil cores. One pump applied the rainfall to the top of the soil core, the second pump sucked the leachate at the bottom of the core. The rainfall intensity, and thus also the flow, were regulated through setting the speed of the pump and the selection of identical 1.3-mm internal diameter tubing.

The system was hermetically closed to maintain the suction. According to the calibrations before and after the experiments rainfall was applied at a rate of 5.3 ± 0.03 mm h-1. Furthermore, the flow was checked by regularly weighing the leachate samples collected automatically in 8-min time intervals.

Before the actual experiments two to three pore volumes of a 0.01 M CaCl2-solution were applied to the soils until steady-state flow conditions and a constant electrical conductivity (determined online) of the leachate samples were reached. We then applied ~10% to 20% of one pore volume (~84 mL) of an aqueous solution containing bromide and Cu at 460 g L-1 KBr and 192 mg L-1 CuN2O6.

The Cu application rate represented twice the recommended maximum annual Cu application rate of 1.2 kg ha-1. The pulse was followed by two to three pore volumes of 0.01 M CaCl2-solution until the electrical conductivity reached its initial value. The conditions of the transport experiments are summarised in Table 2.

Every 10th leachate sample was divided into two sub-samples for analysing Cu and bromide. The leachate samples for Cu analysis were acidified with 0.1% HNO3 (Fisher ultrapure grade 67-69%).

All leachate samples were stored in the dark at 4°C until analysis. 1998). The limit of detection was 4 µg L-1 and the limit of quantification was 7 µg L-1 at a signalto-noise ratio of 3. The other sub-sample of each leachate sample was analysed for Cu at ISTerre laboratory, Grenoble, France via inductively coupled plasma-atomic emission spectrometry (ICP-AES Varian 720-ES Optical Emission Spectrometer) with external calibration. Standards were prepared with a Multi-Element Standards for ICP (Roth) solution. The limit of detection was 0.5 µg L-1 and the limit of quantification was 1 µg L-1 at a signal-to-noise ratio of 3.

Directly after completion of the experiment, the soil cores were weighed to determine the soil's saturated water content. The cores were sliced into 10 segments; from 0-3 cm depth the intervals were 0.6 cm, from 3-6 cm depth the slices were 1 cm, and the remaining soil was sliced into two 2.5-cm slices. The total Cu distribution in the soil with depth for these soil cores, as well as for three additional soil cores from each orchard, were analysed using a standard nitric hydrochloride acid digestion method (US EPA 3050B). In brief, each soil layer was dried at 105°C for 24 h. Then, the soil was manually ground in a ceramic mortar. The Cu determination was completed in triplicate for each soil layer. An aliquot of 1 g of soil was added to 5 mL of aqua regia (HNO3 : 2 HCl) in a 50-mL glass tube. The samples were incubated for 12 h on a temperature adjustable hot plate at 70°C. After diluting the samples with 50 mL of MilliQ water they were filtered with a 0.45 µm cellulose acetate syringe filter (Millipore). A droplet of 1 M HCl was added to preserve the samples until analysis. The analysis was carried out in triplicate for each soil layer by ICP-AES as described above.

In addition, total nitrogen (N) and C concentrations in subsamples of the soil layers were determined by the non-purgeable organic carbon method by 680°C combustion catalytic oxidation in a Shimadzu TOC-V CSN device coupled to a non-dispersive infrared NDIR detector. Standard solutions were prepared from potassium biphthalate and ultrapure water (MilliQ). The pH of the soil layers was also determined. About 10 g of sieved soil was mixed with 10 mL of Milli-Q water and the suspension was stirred for 5 min. The suspension was then left to settle for ~1 h and the pH was measured in the supernatant.

Modelling of transport processes

We used the software 'OM-MADE' to fit a numerical solution of the advection dispersion equation (ADE) coupled with the mobile-immobile model or a mobile-mobile model to the experimental data.

The simulation model was developed for the flow through karst material [START_REF] Collon | Simulation of multiple-peaked skewed breakthrough curves with the new open-source program OM-MADE[END_REF]. It is under publication in the Github repository (https://github.com/OM-MADE/OM-MADE). This is the first application of the software to soils.

The governing equations of the mobile-immobile transport model, which is the classical physical non-equilibrium model assuming two distinct liquid regions in the soil, a mobile region with flowing water and a region with immobile water (van Genuchten and Wagenet 1989), are as follows for nonadsorbing solutes:

  2 * 2 m m m m im m m C Q C C D C C LC t A x x              (5)   im im m im C CC t     (6)
where θ m is the mobile water content, C m is the solute concentration in this mobile water (kg m -3 ), Q is the flux (m 3 s -1 ), A is the total area of the soil core (m 2 ), D* is the dispersion coefficient (m 2 s -1 ), denoting the effective dispersion in the soil, α is the exchange rate between the mobile and the immobile water phases (s -1 ), C im is the solute concentration in the immobile water (kg m -3 ), and L is a loss rate (s -1 ). The total water content of the soil, θ is the sum of the mobile water content, θ m and the immobile water content, θ im . The outflow concentration C out (kg m -3 ) equals the solute concentration of the mobile water fraction, C m .

The mobile-mobile transport model similarly assumes two distinct liquid regions in the soil but both have mobile water phases, which have different flow velocities. Accordingly, the governing equations of the model can be written as:

  2 * 1 1 1 1 1 1 2 1 2 C Q C C D C C t A x x             (7)   2 * 2 2 2 2 2 2 1 2 2 C Q C C D C C t A x x             (8)
where C 1 and C 2 are the solute concentrations (kg m -3 ) in the fast and slow mobile phases, Q 1 and Q 2 are the fluxes (m 3 s -1 ) (with the total flux Q = Q 1 + Q 2 ) with the fast mobile water content θ 1 and the slow mobile water content θ 2 respectively. D 1 * and D 2 * are the dispersion coefficients (m 2 s -1 ) for the fast and slow flow. Here, the total soil water content θ is the sum of the water contents, θ 1 and θ 2 . The outflow concentration C out (kg m -3 ) can be calculated as:

12 12 out QQ C C C QQ  (9)
To select the model that fit the experimental data best we calculated the normalised root mean square error (NRMSE), where lower values indicate less residual variance and thus a better fit:

  2 1 1 n ii i NRMSE s y ny    (10)     1 1 1 1 1 1 2 n ii ii i n yy y x x xx         (11)
where x represents the experimental pore volume and y represents the normalised bromide concentrations measured in the leachate samples, s are the simulated values at the same pore volumes as the measured data, n is the number of experimental measurements, and 𝑦 ̅ is the average experimental bromide concentration, calculated by integrating the experimental curve using the trapeze method.

The 5% arrival time was calculated as an indicator of preferential flow [START_REF] Koestel | Meta-analysis of the effects of soil properties, site factors and experimental conditions on solute transport[END_REF][START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF]). It corresponds to the time between the pulse application and the arrival of the first 5% of the tracer at the bottom of the soil cores normalised to the first temporal moment of the breakthrough curves (BTCs). We also used the holdback factor, which is defined as 'the amount of original fluid remaining in the column when one (water-filled) pore volume of displacing fluid has entered' [START_REF] Koestel | Meta-analysis of the effects of soil properties, site factors and experimental conditions on solute transport[END_REF]). The holdback factor was calculated as the integral of the dimensionless cumulative BTC between zero and one. Absolute experimental data were used for these calculations because flow rate and core size were kept constant in all experiments.

Statistical analysis

The dimensions of the dataset were reduced by principal components analysis (PCA), as implemented in R (R Core Team 2017) using the function prcomp [START_REF] Kassambara | Factoextra: extract and visualize the results of multivariate data analyses[END_REF]. In brief, the orthogonal transformation allowed the visualisation of the potential correlations between the hydrodynamic and soil parameters as well as the macropore characteristics derived by X-ray CT. It also permitted the identification of the minimal number of uncorrelated parameters, the principal components, and ranking them with regard to their contribution to the variance of the data. The principal components were used as new coordinate system, in which each individual core of the dataset was represented using the score value obtained from the PCA. The individual cores in the new coordinate system were then grouped into clusters according to their distances.

In addition, Pearson's correlation analysis and multiple stepwise regression analysis with backward selection were performed with Statistica (StatSoft Inc., Tulsa, USA) to evaluate linear relationships between soil properties, solute transport parameters and soil structural parameters derived from X-ray CT images.

Results and discussion

Effect of orchard management on basic soil properties

The silt loam topsoils in the vine rows under the different orchard management practices had significantly (P < 0.001) different SOC and N contents, which are shown with depth in Figs 1a and1b. Regular compost applications to soils have been shown to result in higher SOC contents in soils [START_REF] Deurer | The impact of soil carbon management on soil macropore structure: a comparison of two apple orchard systems in New Zealand[END_REF][START_REF] Brown | Changes in soil properties and carbon content following compost application: results of on-farm sampling[END_REF][START_REF] Yu | Effects of long-term compost and fertilizer application on stability of aggregate-associated organic carbon in an intensively cultivated sandy loam soil[END_REF]). In our study, the different management of the vine rows, grassed vine rows in the organic orchard versus vegetation-free herbicide-applied rows in the integrated orchard, might also have contributed to the difference observed. Further, regular Cu applications and associated higher Cu concentrations in the O topsoil compared with the I topsoil (Fig. 1d) might have inhibited microbial activity, and thereby reduced C mineralisation rates.

The soil pH was similar in both topsoils (Fig. 1c). The total porosity of the O topsoil was 0.716 cm3 cm-3, significantly (P = 0.01) larger than that of the I topsoil (0.62 cm3 cm-3) (Table 2). The C content was strongly and positively correlated with total porosity (r = 0.95; P = 0.004).

Effect of orchard management on pore architecture

Macroporosity tended to be higher in I (an average of 0.064 cm3 cm-3) than O (0.043 cm3 cm-3) but this difference was not significant (P = 0.2). The macroporosity is within the range of macroporosities reported in the past for arable, horticultural and pastoral soils with similar X-ray CT resolutions [START_REF] Peth | Three-dimensional quantification of intra-aggregate pore-space features using synchrotron-radiation-based microtomography[END_REF][START_REF] Deurer | The impact of soil carbon management on soil macropore structure: a comparison of two apple orchard systems in New Zealand[END_REF][START_REF] Paradelo | X-ray CT-derived soil characteristics explain varying air, water, and solute transport properties across a loamy field[END_REF][START_REF] Müller | Characterising and linking X-ray CT derived macroporosity parameters to infiltration in soils with contrasting structures[END_REF]. The macropore networks are visualised in Fig. 2. We found a weak negative correlation between SOC and macroporosity (r = -0.39; P < 0.05). In contrast, [START_REF] Deurer | The impact of soil carbon management on soil macropore structure: a comparison of two apple orchard systems in New Zealand[END_REF] reported a significant positive correlation between SOC and macroporosity for topsoils of apple orchards under organic and integrated orchard management. It is worthwhile to mention that we found a strong positive correlation between SOC and sub-X-ray CT resolution porosity as mentioned above. [START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF] and [START_REF] Paradelo | X-ray CT-derived soil characteristics explain varying air, water, and solute transport properties across a loamy field[END_REF] found no correlation between CT-derived macroporosity and SOC contents. These inconsistent findings might be explained by different soil types, clay contents and mineralogy, and the role of clay on soil structure [START_REF] Naveed | Revealing soil structure and functional macroporosity along a clay gradient using X-ray computed tomography[END_REF].

The I topsoil had with 43 407 (±6207) macropores per core nearly significant (P = 0.06) more macropores than the O topsoil with 26 988 (±6179) macropores per core. This resulted in macropore densities ranging between 0.07 and 0.1 mm-3 for I and between 0.03 and 0.07 mm-3 for O. These are lower than the macropore density reported for a pastoral Andosol topsoil [START_REF] Müller | Characterising and linking X-ray CT derived macroporosity parameters to infiltration in soils with contrasting structures[END_REF]) but an order of magnitude higher than those reported by [START_REF] Zhang | Quantification of 3D macropore networks in forest soils in Touzhai valley (Yunnan, China) using X-ray computed tomography and image analysis[END_REF] for forest soils. The average macropore diameter was significantly (P < 0.05) larger for the O topsoil (0.64 ± 0.0032 mm) than for the I topsoil (0.59 ± 0.0018 mm). Assuming that a large proportion of macropores were formed through root activities [START_REF] Hatano | Relationship between the distribution of soil macropores and root elongation[END_REF], our results of more and smaller macropores in I compared with O are consistent with previous research on the effect of orchard groundcover management on tree roots. The longevity of shallow tree roots in orchard rows kept weed-free through herbicide applications has been reported to be lower than in grassed orchard rows, leading to smaller root diameters [START_REF] Yao | Apple root growth, turnover, and distribution under different orchard groundcover management systems[END_REF]). This might explain the smaller average macropore diameter found in I compared with O. In addition, topsoil temperatures are higher in tree rows that are kept weed-free than in grassed rows [START_REF] Merwin | Orchard groundcover management impacts on soil physical properties[END_REF], which might affect root growth. Kiwifruit vine root growth is limited above 25°C [START_REF] Smith | Effect of oxygen supply and temperature at the root on the physiology of kiwifruit vine[END_REF], which might also have contributed to the smaller average macropore diameter in I compared with O. In addition, consistent with our findings of more macropores in I than in O, [START_REF] Yao | Apple root growth, turnover, and distribution under different orchard groundcover management systems[END_REF] observed more shallow roots in weed-free apple orchard tree rows than in grassed tree rows. This was explained by lower competition for water and nutrients in weed-free rows than in grassed rows [START_REF] Parker | Orchard floor management affects peach rooting[END_REF]). Moreover, wetting-drying cycles also contribute to macropore formation [START_REF] Horn | Structure formation and its consequences for gas and water transport in unsaturated arable and forest soils[END_REF]. Wetting-drying cycles would have occurred more often in the weed-free vine rows of I than O due to larger temperature and moisture gradients [START_REF] Ananyeva | Can intra-aggregate pore structures affect the aggregate's effectiveness in protecting carbon?[END_REF]. Furthermore, macro-fauna and earthworms in particular are important soil engineers. For example, in a catchment study earthworm numbers were found to be positively correlated with the number of macropores [START_REF] Van Schaik | Linking spatial earthworm distribution to macropore numbers and hydrological effectiveness[END_REF]. Any of these reasons could explain our findings of more and on average smaller macropores in I than in O.

Average tortuosity was 2.37 ± 0.1 in O and compared with 2.47 ± 0.03 in I (P = 0.18; Fig. 3a). It was comparable to the tortuosity of an Andosol topsoil under pasture [START_REF] Müller | Characterising and linking X-ray CT derived macroporosity parameters to infiltration in soils with contrasting structures[END_REF] with values higher than those reported for arable soils [START_REF] Katuwal | Effects of CT number derived matrix density on preferential flow and transport in a macroporous agricultural soil[END_REF]. The connectivity density was significantly (P < 0.05) lower in O (average of -0.23 ± 0.43) than in I (0.19 ± 0.06). The lower the connectivity density value, the higher the pore connectivity (Fig. 3b). Although the total number of pores was significantly (P < 0.05) higher in O than in I (Fig. 3c), the number of connected pores in the two orchard soils was about the same (P = 1). The pore shapes were similar in both soils, and were dominated by tubular pores, which are indicative of pores formed by roots, supporting the assumption that roots were important for macropore formation in these orchards. Macropores in I tended to be more elongated (Fig. 3d).

Effect of orchard management on hydrodynamic soil properties

The dimensionless measured and simulated BTCs for the inert tracer bromide through the two orchard topsoils are shown in Fig. 4. The pulse volume between the six cores varied between 83.2 and 87.4 mL with an average of 84.3 ± 0.64 mL, which can be explained by the difficulties of establishing the same fluxes in the soil cores. The experimental parameters are summarised in Table 2. Table 3 presents the corresponding modelling parameters. The Peclet number for all experiments was >5, which according to [START_REF] Kutílek | Soil hydrology[END_REF] indicates that molecular diffusion processes can be neglected and the hydrodynamic dispersion can be assumed to be equal to the advective dispersion. Dispersivity (λ; m) can then be estimated as D*/v with v being flow velocity (m s-1). Dispersivity varied between 3 × 10-3 and 8 × 10-3 m for the O topsoil and was in the same range as reported in the literature for intact Andosol soil samples (e.g. [START_REF] Magesan | Solute movement through an allophanic soil[END_REF][START_REF] Prado | Transport, sorption and degradation of atrazine in two clay soils from Mexico: Andosol and Vertisol[END_REF].

Dispersivity displayed a much wider range for the I topsoils (Table 3).

All BTCs from the O topsoil cores were more or less symmetrical with the centre of gravity at or close to one pore volume (Figs 4a-c). Bromide sorption was excluded based on complementary batch experiments (data not shown). The ADE-model fit the measured data reasonably well for all three soil cores, with NRMSE of 5.6%, 4.6% and 5.3% for O1, O2 and O3 respectively. For O1, this NRMSE was achieved by adjusting the bromide input to address the poor mass balance, which might be related to an experimental problem with the bromide analysis. In spite of this issue, the results suggest that flow processes through the O topsoils were dominated by homogeneous flow through a wellconnected pore network. In contrast, the BTCs from I were asymmetrical and showed a tailing characteristic of preferential flow [START_REF] Brusseau | Modeling solute transport in structured soils: a review[END_REF]. The bromide peak occurred in all three cores well before one pore volume ( Figs 4d-f). The experimental data fitted well with the mobileimmobile model for two of the soil cores collected from the integrated orchard (I1 and I3) with NRMSE-values of 2.7% and 5.2%. For I2, the mobile-immobile model resulted in an NRMSE of 8.2%. The BTC for I2 was characterised by a double peak with the first larger peak before one pore volume. This transport behaviour was best simulated with the mobile-mobile model (NRMSE, 3.8%). This model implies two mobile water phases, which in our case had similar velocity but differed in dispersivity (Table 3). This difference in dispersivity was the main cause for the double-peak behaviour. Applying the mobile-immobile or mobile-mobile models greatly improved the simulations of the BTCs for I1 and I3. The respective NRMSE-values using the ADE model were 18%, 25% and 6% respectively. The modelled hydrodynamic parameters for the different modelling approaches are summarised in Table 3. The mobile water contents for I1 and I3 were 0.01 and 0.05, and thus very low, indicating that most of the water was immobile. Usually, values >0.5 are reported. However, such low mobile water fractions are not unrealistic but highlight the importance of macropore flow [START_REF] Griffioen | Interpretation of two-region model parameters[END_REF]. The first order mass transfer coefficient between mobile and immobile zone, α was 0.89 h-1 for I1 and 1.38 h-1 for I3, suggesting that the exchange in I3 was quicker than in I1.

The 5% arrival time of bromide was positively correlated with SOC content (r = 0.83; P < 0.05), while the holdback factor was negatively correlated with SOC (r = -0.95; P < 0.05). The 5% arrival times for the O topsoils were with 0.34, 0.37 and 0.39 for O1, O2 and O3 respectively, and thus close to the suggested threshold for advective-dispersive transport of 0.4 [START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF] supporting the findings of our OM_MADE simulations using the ADE equation. [START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF] reported similar relationships of less preferential flow in soils with higher C contents for soils from a single field with a SOC gradient from 40 to 160 g kg-1. The respective 5% arrival times for the I topsoils were 0.113, 0.153 and 0.283 for I1, I2 and I3 respectively and significantly shorter (P = 0.03) than those for O.

Relationships between X-ray CT parameters and flow

Nine of the generic soil parameters and 14 hydrodynamic parameters derived from the transport experiments with bromide were considered to be relevant for flow processes through soils (Table 4).

To discern the relationships between structural and functional soil parameters and to reduce the number of multiple dependent variables addressing the problem of multicollinearity, we conducted a PCA analysis. For the PCA, we considered structural parameters for the percolating pores only because only pores well connected in the vertical direction are considered relevant for preferential flow [START_REF] Jarvis | Connectivity and percolation of structural pore networks in a cultivated silt loam soil quantified by X-ray tomography[END_REF]. Four principal components had eigenvalues >1. The first two components with eigenvalues >1 were extracted, which explained ~75.4% of the total variance of the original data (Table 4, Fig. 5). The first component explained the highest variance in the results (54.5%), whereas the second principal component accounted for 20.9% of the variance. The relation with the parameters is visualised as vectors in Fig. 5a, and the position of each core is shown in the new 2D-coordinate system in Fig. 5b. The vector loading plot (Fig. 5a) shows the parameters with the greatest effect on discriminating the orchards. The PCA separated the two orchards along PC1 (Fig. 5b). The soils from O clustered at the left side of the biplot, and were clearly separated from the soils from I, which were further spread apart from each other on the right side of the biplot. The soils from O were separated from the I topsoils along PC1 due to larger total porosity, mobile water content, pore diameter, SOC content, and lower holdback factor and bulk density (Fig. 5a). For I, a spread along PC2 was observed, which was related to tortuosity, pore curvature, sphericity and the pore water velocity. [START_REF] Koestel | Evaluation of nonparametric shape measures for solute breakthrough curves[END_REF] conducted a meta-analysis with 733 BTCs in order to improve the understanding of relationships between soil properties and site factors and the occurrence of preferential flow. They identified the 5% arrival time and the holdback factor as useful indicators for the strength of preferential flow. We conducted multiple linear regression analyses to identify the models that explained 5% bromide arrival time and its holdback factor best including both CT-derived structural and general soil parameters. The best and most significant model (P = 0.001; Fig. 6a), which explained ~98% of the variability of the 5% bromide arrival time included SOC (P = 0.0005) and tortuosity (P = 0.002). The parameters C (P < 0.001), connectivity (P = 0.006) and tortuosity (P = 0.009) were selected in the best model for the holdback factor (P < 0.001; Fig. 6b). These results highlight the importance of SOC contents for the occurrence of preferential flow, as also recently highlighted by [START_REF] Larsbo | Preferential transport in macropores is reduced by soil organic carbon[END_REF]. They also support previous findings that structural soil parameters can considerably improve predictions of preferential flow [START_REF] Luo | Quantitative relationships between soil macropore characteristics and preferential flow and transport[END_REF][START_REF] Katuwal | Effects of CT number derived matrix density on preferential flow and transport in a macroporous agricultural soil[END_REF][START_REF] Paradelo | X-ray CT-derived soil characteristics explain varying air, water, and solute transport properties across a loamy field[END_REF].

Effects of orchard management on Cu leaching and accumulation

We did not detect Cu in the leachate samples of O. For I, we detected traces of Cu (concentrations close to the detection limit of 0.5 µg L-1) but only in a few leachate samples collected before one pore volume. Copper was probably preferentially transported in I, albeit at very low concentrations.

Quantification was impossible with our analytical protocol.

We analysed the total Cu distribution with depth in the two orchard soils as well as in the cores after completion of the transport experiments. The dotted lines in Fig. 7a present the average total Cu concentration in the topsoils before the transport experiments. The concentrations were significantly (P > 0.05) higher in O than I. This was explained by the fact that orchardists following organic management principles traditionally rely on Cu as fungicide. The Cu enrichment factor for the transport experiment is presented in Fig. 7b; the applied Cu accumulated in the top 2.5 cm in O, whereas it was moved at least to 8 cm depth in I. This is supported by the finding of Cu traces in the leachate samples of I. About 85% of the variability of the total Cu concentrations in the segments of the two topsoils was explained by SOC, pH, pore circularity and equivalent pore diameter (Fig. 8; adjusted R2 = 0.85).

Conclusions

Our study highlights the significant effect of management on SOC content and network architecture of pores >92 µm in topsoils under integrated and organic kiwifruit production. As expected, SOC contents were significantly higher in O than in I topsoils. Macropores were larger with a higher connectivity in the O topsoil than in the I topsoil, but they were significantly fewer in number.

Predominantly tube-shaped macropores identified root activities as important factor for macropore formation. In accordance with this, observed differences in the macropore network may be attributed to the effect of herbicides on shallow kiwifruit roots. Significant relationships between soil parameters and indicators of the strength of preferential flow of our water tracer (bromide) verified that structural parameters derived by X-ray CT might be useful for predicting the functional solute transport parameters 5% arrival time and the holdback factor. In addition, predicting the depth distribution of the strong-sorbing solute Cu was best when combining SOC contents and pH with the structural parameters macropore circularity and diameter. Organic orchard management improved the filtering function of the topsoil through increasing SOC content and changing pore architecture. It has to be noted that our results are based on a very limited number of soil cores from only two orchards and need to be confirmed by analysing a larger number of samples.

  Bromide concentrations were quantified by ion chromatography (Metrohm 732/733 separation centre, Metrosep A Supp 16-150 mm, and Metrosep C2-150 mm anion and cation columns) based on the USEPA Method 300.0 (Pfaff

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 12345678 Fig. 1. (a) Total carbon concentration, (b) total nitrogen concentration, (c) pH and (d) total copper concentration with depth determined for the soil cores from the organic and the integrated kiwifruit orchards. Presented are mean values and standard error of the means for depth intervals and three cores per orchard.

Table 1 . Generic soil profile description of the organic and integrated kiwifruit orchards located in the Bay of Plenty, New Zealand

 1

	Orchard	Depth (m)	Texture	Sand (%) Silt (%)	Clay (%) Carbon	Bulk	density
							(%)	(10 3 kg m -3 )
	Organic	0-0.3	Silt loam	41.34	35.64	23.02	5.94	0.865
	orchard	0.3-0.5		42.84	41.64	15.52	2.87	0.936
		0.5-1		54.32	33.98	11.70	1.83	0.943
	Integrated	0-0.3	Silt loam	46.40	31.64	21.96	2.46	0.987
	orchard	0.3-0.5		43.40	40.14	16.46	1.41	0.961
		0.5-1		55.90	29.14	14.96	0.51	1.014

Table 2 . Parameters of the transport experiments with copper and the water tracer bromide through 111 mm-long intact soil cores with an ID of 104 mm collected from the organic (O) and the integrated (I) kiwifruit orchard in the Bay of Plenty, New Zealand

 2

	Core	Bulk	Total	Pore		Vol. water	Pore	water	Pulse duration in pore
	ID	density (m 3	porosity (m 3	volume		content (m 3	velocity (m s -1 )	volume
		m -3 )	m -3 )	(m 3 )		m -3 )
	O1	0.81	0.696	5.77	×	0.624	3.04 × 10 6	0.144
				10 4
	O2	0.7	0.745	5.24	×	0.56	2.55 × 10 6	0.160
				10 4
	O3	0.78	0.706	5.45	×	0.586	2.94 × 10 6	0.154
				10 4
	I1	1.02	0.615	4.37	×	0.472	3.60 × 10 6	0.191
				10 4
	I2	0.99	0.626	4.85	×	0.509	3.42 × 10 6	0.180
				10 4
	I3	1.006	0.620	4.70	×	0.506	3.67 × 10 6	0.178
				10 4

Table 3 . OM-MADE modelled hydrodynamic parameters, variables for the transport experiments with bromide in the organic (O) and the integrated (I) kiwifruit orchard in the Bay of Plenty, New Zealand. The following model types were applied: (1) advection-dispersion equation (2), mobile-immobile model, and (3) mobile-mobile model. Results with the lowest normalised root mean square error (NRMSE) are shown

 3 In the mobile-mobile model, θ m is the fast mobile water content, θ im is the slow mobile water content, PWV is the flow rate in the fast mobile water phase, PWV 2 in the slow mobile water phase, D* is the dispersion coefficient in the fast mobile water phase, and D 2 * is the dispersion coefficient in the slow

	mobile water phase
	Core	Mobile	Immobile	Pore	Pore	Dispersion	Dispersion	Dispersivity	Exchange	NMRSE	Model
	ID	water content θ m (m 3 m -3 )	water content θ im (m 3 m -3 )	velocity PWV (m s -1 )	velocity PWV 2 (m s -1 )	coefficient , D* (m 2 s -1 )	coefficient s -1 ) , D 2 * (m 2	, λ (10 3 m)	(s -1 ) rate, α	(%)
	O1	0.624	0	3.1 × 10 6	-	2.6 × 10 8	-	8.38	-	5.6	1
	O2	0.560	0	2.55 × 10 6	-	7.47 × 10 9	-	2.93	-	4.6	1
	O3	0.586	0	2.94 × 10 6	-	1.46 × 10 8	-	4.99	-	5.3	1
	I1	0.012	0.460	1.39 × 10 4	-	1.83 × 10 8	-	0.13	1.15 × 10 2	2.7	2
	I2	0.215	0.294	4.23 × 10 6	3.09 ×	7.15 × 10 8	4.9 × 10 9	16.9 (1.58)	0	3.8	3
					10 6
	I3	0.052	0.454	3.54 × 10 5	-	3.06 × 10 8	-	0.86	1.75 × 10 2	5.2	2

Table 4 . Correlations between principal components PC1, PC2, PC3 and PC4 and the parameters included in the PCA. Boldfaced values are >0.7

 4

	Factor 1	Factor 2	Factor 3	Factor 4

Acknowledgements

This work was started during a secondment of the first author at the Institut des Géosciences de l'Environnement (IGE) in Grenoble as invited scientist by Observatoire des Sciences de l'Univers de Grenoble (OSUG). The continuation of the work was supported by the New Zealand-European Union (EU) Marie Sklodowska-Curie Research and Innovation Staff Exchange (RISE) Counterpart Fund from Government funding (Contract 645717), administered by the Royal Society of New Zealand.

Article published in Soil Research (www.publish.csiro.au/SR/SR18293) doi.org/10.1071/SR18293

Conflicts of interest

The authors declare no conflicts of interest.