Marco Dentz
email: marco.dentz@csic.es

Alessandro Comolli

Vivien Hakoun

Juan Hidalgo

Transport upscaling in highly heterogeneous aquifers and the prediction of tracer dispersion at the Macrodispersion Experiment (MADE) site

We present an upscaled Lagrangian approach to predict the plume evolution in highly heterogeneous aquifers. The model is parameterized by transport independent characteristics such as the statistics of hydraulic conductivity and the Eulerian flow speed. It can be conditioned on the tracer properties, and flow data at the injection region. Thus, the model is transferable to different solutes and hydraulic conditions. It captures the large scale non-Gaussian features for the evolution of the longitudinal mass distribution observed for the bromide and tritium tracer plumes at the macrodispersion experiment (MADE) site (Columbus, Mississippi, USA), which are characterized by a slow moving peak and pronounced forward tailing. These large scale features are explained by advective tracer propagation due to a broad distribution of spatially persistent Eulerian flow speeds as a result of spatial variability in hydraulic conductivity.

I. INTRODUCTION

The upscaling and prediction of tracer transport in highly heterogeneous porous and fractured media is of central concern in a broad range of subsurface applications from groundwater management to underground gas and waste storage. This task is challenging due to strong spatial variability in hydraulic conductivity values encountered in geological media [START_REF] Bear | Dynamics of Fuids in Porous Media[END_REF].

Spatial variability in hydraulic conductivity induces spatial fluctuations in the flow and transport velocities, whose impact on large scale tracer migration has been quantified in terms of macrodispersion coefficients [START_REF] Gelhar | [END_REF]3]. For moderately heterogeneous media, Gelhar and Axness [START_REF] Gelhar | [END_REF] used a stochastic modeling approach [START_REF] Dagan | Flow and transport in porous formations[END_REF][START_REF] Gelhar | Stochastic subsurface hydrology[END_REF][START_REF] Rubin | Applied stochastic hydrogeology[END_REF] to express the longitudinal macrodispersion coefficients in terms of the mean flow velocity, and the correlation length and variance of the logarithm of hydraulic conductivity. This expression is a central results for hydrogeological prediction because it allows to forecast macroscopic transport features based on transport independent observables. However, spatial heterogeneity gives rise to transport behaviors that can be very different from the ones predicted by advection-dispersion models characterized by constant macrodispersion coefficients. For pointlike solute injections, the latter predicts Gaussian shaped tracer plumes and breakthrough curves, while observed distributions are typically found to be non-Gaussian [START_REF] Adams | [END_REF][8][9][10][11]. This is the case for the tracer plumes monitored during the macrodispersion experiments conducted in the alluvial aquifer underlying the Columbus Air Force Base in northeastern Mississippi [START_REF] Adams | [END_REF]12]. Spatial tracer distributions are characterized by strongly non-Gaussian shapes characterized by a slowly moving peak and a pronounced forward tail. These and other observations of anomalous solute dispersion have spurred the development of non-Fickian transport theories [13][14][15][16]. Those include multirate mass transfer approaches [17,18], continuous time random walks [13,19], fractional advection-dispersion models [20,21], time-domain random walks [22,23], and space-time non-local advection-dispersion equations [24,25]. While all of these approaches provide dynamic frameworks to model non-Gaussian large scale transport features, their parameterization in terms of transport independent parameters and thus their predictive power remain open questions.

The migration of the tracer plume at the MADE site was modeled by Berkowitz and Scher [26] with the continuous time random walk (CTRW) approach based on an empirical distribution of transition length and times, which reflects a broad distribution of mass transfer time scales. Benson et al. [27], Schumer et al. [28] and Zhang and Benson [29] modeled the tracer plumes of the MADE-1 experiment using a fractal mobile-immobile model that accounts for both solute retention due to a broad distribution of mass transfer time scale, and preferential transport due to a broad distribution of mass transfer length scales. The characteristic exponents of the model are estimated from the experimental data. Harvey and Gorelick [30] modeled the plume evolution using a multirate mass transfer (MRMT) model that considers rate limited mass transfer between a mobile domain with accelerating flow, and an immobile zone, which represents intragranular porosity, low permeability zones, dead end pores and surface sorption sites. These different processes are represented by a capacity coefficient, a rate coefficient, a retardation coefficient, a velocity parameter and an acceleration parameter, which are estimated from the experimental mass data. These modeling approaches propose a range of mass transfer processes and invoke broad distributions of mass transfer time scales and lengths in order to simulate the impact of medium and flow heterogeneity on large scale transport. A key question in order to constrain such models refers to the dominant local scale mechanisms that cause non-Fickian large scale transport.

Mobile-immobile mass transfer has been studied as a mechanism to explain the decay of the apparent solute mass observed at the MADE-1 and MADE-2 experiments [30]. Feehley et al. [31] and Guan et al. [32] combine a stochastic reconstruction of aquifer heterogeneity based on fractal Brownian motion with mobile-immobile mass transfer in order to assess the plume evolution in the MADE-2 experiment. A series of works [33][34][35] have shown that a local scale advection-dispersion model based on the detailed knowledge of the spatial distribution of hydraulic conductivity allows reproducing the spatial distribution of tracer plumes observed at the MADE site. This implies that some key features of transport may indeed be understood from the spatial variability of hydraulic conductivity and thus advective heterogeneity. The question arising from this observation is whether characteristics of non-Fickian transport in strongly heterogeneous media can be predicted based only on a few geostatistical medium characteristics similar to the prediction of macrodispersion in moderately heterogeneous media. Along these lines, several authors [23,[36][37][38] investigated the relation between advective travel times and hydraulic conductivity, and the impact of broad distributions of the logarithm of hydraulic conductivity on non-Fickian transport. For example, the stochastic multi-indicator approach of Fiori et al. [23] was adapted in Fiori et al. [39] to predict the plume evolution of the MADE-1 experiment based on analytical expressions for the advective travel times over rectangular inclusions.

In this paper, we adress the questions of upscaling and prediction of dispersion in highly heterogeneous aquifers in terms of medium and flow heterogeneity, using transport independent parameters. We use a stochastic time domain random walk approach [40,41] that accounts for heterogeneous advection. This approach propagates particles in space and time according to a velocity Markov model that is determined by the distribution of the Eulerian flow speed and correlation length. We discuss the model and its parameterization in terms of the medium and flow properties, and use it for the modeling and interpretation of the one-dimensional tracer profiles of the MADE-1 experiment.

II. UPSCALED TRANSPORT MODEL

We employ the stochastic time-domain random walk (TDRW) implementation presented in [41] in order to build an upscaled model to predict the concentration evolution in highly heterogeneous aquifers, based on the statistical characteristics of the logarithm of hydraulic conductivity, the flow speed, porosity and the mean hydraulic gradient. The method quantifies particle motion in spatially variable flow fields through a stochastic model for equidistantly sampled particle speeds along streamlines. Thus, it propagates solute particles at constant space increments whose duration is obtained from a Markov chain representation for the particle speed. The model is predictive because it can be parameterized by transport independent medium and flow characteristics. This approach has been used for the upscaling and prediction of solute transport from the pore to Darcy scale, and from the Darcy to the field scale in two-and three-dimensional porous and fractured media [11,[41][42][43][44][45][46][47][48]. In the following, we detail the model, its assumptions, and parameterization.

A. Stochastic advective motion

We choose a Lagrangian framework [START_REF] Dagan | Flow and transport in porous formations[END_REF] and focus on the advective motion of tracer particles, which constitute the tracer plume. The temporal variation patters of particle velocities are in general rather complex. They are intermittent because their evolution is governed by a broad distribution of times scales, which are determined by the flow speeds and characteristic length scales imprinted in the medium structure [40,49]. The complexity of intermittent temporal series of particle velocities can be removed by adopting an equidistant point of view. Therefore, unlike classical Lagrangian approaches we do not consider particle motion isochronically, but equidistantly along streamlines. This means, we quantify motion through transitions over constant space increments instead of constant time increments. It has been shown for the pore and the Darcy scale [40,42,48,50] that equidistant series of particle velocities can be modeled as Markov chains.

Thus, advective particle motion is quantified by the following stochastic evolution equations for the longitudinal particle position x(s) and the advective travel time t a (s) [41,49],

dx(s) ds = χ -1 , dt a (s) ds = 1 v (s) ( 1
)
where s is the distance along streamline, and v (s) > 0 the particle speed, this means the magnitude of the equidistant particle velocity. The distance s along the streamline is in general larger than the linear travel distance x in the mean flow direction. They are projected onto each other through the advective tortuosity χ, which is given by the ratio of the average flow speed q and the average velocity component q 1 in mean flow direction, χ = q / q 1 [41].

The angular brackets denote suitably chosen averages, here spatial and ensemble average. The particle speed v (s) describes an ergodic and stationary Markov process characterized by the steady state distribution p v (v). The latter is a Lagrangian quantity, which, however, can be related to the distribution of Eulerian flow speed p e (v) via the flux-weighting relation [41,49]

p v (v) = vp e (v) v , (2)
Note that v (s) evolves as a function of distance s along a streamline, while in classical Lagrangian methods [START_REF] Dagan | Flow and transport in porous formations[END_REF], particle speeds evolve with travel time. This difference explains the flux-weighting relationship (2) as discussed in detail in Dentz et al. [49] and Hakoun et al. [40].

The Eulerian velocity v(x) is given in terms of the Darcy flux q(x), porosity φ, and the retardation coefficient θ, which accounts for linear equilibrium sorption. We assume that both porosity and retardation are constant. The Darcy flux satisfies

q(x) = -K(x)∇h(x), ∇ • q(x) = 0, (3)
where K(x) is hydraulic conductivity and h(x) hydraulic head. The flow problem ( 3) is solved numerically for a unit mean hydraulic gradient aligned with the coordinate axis in one-direction, and a unit geometric mean conductivity K g = 1, see Appendix A. The Eulerian velocity v(x) is given in terms of the resulting Darcy velocity q(x) as

v(x) = v 0 q(x) v 0 = K g J/φθ ( 4
)
where J is the magnitude of the mean hydraulic gradient. The velocity v 0 and the correlation length c define the characteristic time τ v = c /v 0 . Note that changes in v 0 affect transport only through a rescaling of time. The Eulerian speed is v(x) = |v(x)|, and its distribution is, as introduced above, denoted by p e (v). The magnitude of q(x) is denoted by q(x). Based on (4), the Eulerian speed distribution is given in terms of the distribution p q (q) of q(x) as

p e (v) = p q (v/v 0 )/v 0 . (5)
In this work, hydraulic conductivity is represented by a three-dimensional lognormally distributed spatial random field. The horizontal correlation length is h = 5 v and the variance of the logarithm f (x) = ln[K(x)] of hydraulic conductivity is set to σ 2 f = 5.9 reflecting the conductivity statistics estimated for the MADE site [51]. We find q 1 = 3.04 and q = 3.5. This implies that the effective conductivity is K e = 3.04K g , and the advective tortuosity χ = 1.15. Figure 1 shows the distributions of f (x) = ln[K(x)] and ν(x) = ln[q(x)]. The distribution p ν (ν) of ν(x) can be modeled by the skewed Gaussian distribution

p ν (ν) = exp -(ν-m) 2 2σ 2 ν 2πσ 2 ν Φ[α(ν -m)/ σ 2 ν ], (6)
where Φ(x) is the cumulative unit Gaussian distribution. The parameters in Figure 1 are σ 2 ν = σ 2 f = 5.9, α = 1.29 and m = 1.47.

Regarding the propagation of the particle speeds v (s) from their initial values v 0 = v (0), the key feature is the spatial persistence of speeds, which is characteristic for transport in steady flow fields [19,22,42,49]. This spatial persistence is at the origin of intermittent Lagrangian flow properties in porous media found both at the pore and Darcy scales [40,45,50,52]. Thus velocity transitions have been modeled, for example through independent sampling after fixed persistence lengths [13], Bernoulli processes [11,40,49], and by empirical velocity transition matrices [42,43]. In this work, we consider the evolution of the normal score w(s) of v (s), which is given by

w = Φ -1 [P v (v)], v = P -1 v [Φ(w)]. (7)
The cumulative probability of p v (v) is denoted by P v (v), its inverse by P -1 v (v). The Doob theorem states that the only process that is both Markovian and Gaussian is the Orrnstein-Uhlenbeck process [53,54]. For this reason, we model the evolution of w(s) through the Ornstein-Uhlenbeck process [54] where η(s) is a Gaussian white noise characterized by 0 mean and correlation η(s)η(s ) = δ(s -s ). The steady state distribution of w(s) is the unit Gaussian. The initial values w(0) of the normal scores are obtained from v(0) through the Smirnov transform [START_REF] Adams | [END_REF]. The initial velocity distribution is denoted by p 0 (v). The correlation length c of the normal score w(s) is set equal to c = 8 h /3 [22,40,41]. This approach has been used and verified for the prediction of equidistant Lagrangian velocity series and intermittent transport characteristics in pore and Darcy scale heterogeneous porous media [40,41,45,50].

dw(s) ds = --1 c w(s) + 2 -1 c η(s), (8)
Equations (1) together with the evolution equation for v (s), describe the propagation of the particle position x(s) and velocity v (s) from the initial values x(0), and v(0), which are distributed according to ρ(x), and p 0 (v). The initial value of time is t(0) = 0. The projected streamwise concentration profile c(x, t) of the mobile, this means non-adsorbed solute in this framework is obtained by

c(x, t) = θ -1 δ [x -s(t)/χ] , (9)
where s(t) = max[s|t a (s) ≤ t]. The angular brackets denote here the average over all particles. The TDRW approach is intrinsically mass conservative. Thus its density describes the total, this means adsorbed and non-adsorbed concentrations, which is the reason why the non-adsorbed concentration is obtained from the total density through division by θ as given in (9). This stochastic TDRW approach is parameterized by the Eulerian speed distribution, advective tortuosity and the correlation length c . These quantities are constrained by the geometric mean conductivity, the variance and correlation lengths of the logarithm of hydraulic conductivity, the mean hydraulic gradient, and porosity, without non-physical fitting parameters. This means the model is parameterized in terms of geostatistical and hydraulic characteristics that can be determined by transport independent measurements, and is in this sense predictive. In addition, the approach can account for linear equilibrium sorption.

In the following, we use this upscaled transport model to predict the longitudinal plume evolution of the MADE-1 and MADE-2 experiments. Its numerical implementation is detailed in Appendix B.

III. PREDICTION OF DISPERSION AT THE MADE SITE

We use the upscaled transport model presented in the previous section for the prediction of the streamwise tracer profiles of the MADE-1 and MADE-2 experiments [START_REF] Adams | [END_REF][START_REF] Boggs | Transport of tritium and four organic compounds during a natural-gradient experiment (MADE-2)[END_REF]. For the MADE-1 experiment, we consider six snapshots of the longitudinal bromide concentration profiles [START_REF] Adams | [END_REF]39]. For the MADE-2 experiment, we consider two snapshots of the longitudinal tritium profiles [34,[START_REF] Boggs | Transport of tritium and four organic compounds during a natural-gradient experiment (MADE-2)[END_REF]. In the following, we first discuss the model parameters for the MADE site, and then report on the prediction of the longitudinal concentration profiles for the MADE-1 and MADE-2 experiments. [START_REF] Boggs | [END_REF] and Harvey and Gorelick [30].

Parameter σ 2 f Kg [m/s] h [m] v [m] J [-] φ θ [Br] θ [ 3 H
A. Hydraulic and transport parameters for the MADE site

We consider purely advective transport, which is a reasonable assumption given that the Péclet number at the MADE site is around 10 3 [39]. In order to parameterize the model, we rely on the description of the experimental conditions in Boggs et al. [12] and Adams and Gelhar [START_REF] Adams | [END_REF] and the geostatistical characterization of the hydraulic conductivity field given in Bohling et al. [51]. Thus, the medium porosity, and the magnitude of the mean hydraulic gradient are set equal to φ = 0.31, J = 3.6 × 10 -3 [START_REF] Adams | [END_REF]12,39]. Bohling et al. [51] estimate for the variance of the logarithm of conductivity σ 2 f = 5.9 and for the horizontal and vertical correlation lengths h = 9.1 m and v = 1.8 m. For the geometric mean conductivity these authors obtain the average value of K g = 6.7 × 10 -6 m/s. The upper and lower limits of the 95% confidence interval around this value are K g = 10 5 m/s and K g = 4.3 × 10 -6 m/s. Here we use the value K g = 5.5 × 10 -6 m/s because it provides a better prediction of the observed peak concentration values. The model predictions for K g = 6.7 × 10 -6 m/s and a discussion of the evolution of the peak values are given in Appendix D.

The values for the predictions of the longitudinal concentration profiles shown below are summarized in Table I. The characteristic velocity v 0 is given by v 0 = 6.39 × 10 -8 θ -1 m/s and the average transport velocity by v 1 = K e J/φθ = 1.94 × 10 -7 θ -1 m/s. We use for the prediction of both the MADE-1 and MADE-2 data the same medium and hydraulic parameters. Boggs and Adams [START_REF] Boggs | [END_REF] found from column experiments the values θ = 1.2 and θ = 1 for the retardation coefficients of bromide and tritium. The retardation coefficient accounts for the fact that only the mobile, non-adsorbed solute is measured [30,[START_REF] Boggs | [END_REF], and in this sense accounts for a measurement bias.

The parameter values specified above determine the propagator of the upscaled transport model. In order to predict the plume evolution, we need both the propagator, and the initial conditions, this means here, the initial particle positions, and initial particle speed distribution. As described in Boggs et al. [12] and Stauffer et al. [START_REF] Stauffer | Degradation of aromatic hydrocarbons in an aquifer during a field experiment demonstrating the feasibility of remediation by natural attenuation[END_REF], in both experiments tracer was injected into five 5.2-cm-diameter injection wells separated by 1 m in a linear array perpendicular to the mean flow direction. In both experiments, tracer was injected in about 10 m 3 of groundwater for the duration of 48.5 hours. Thus, while the initial distribution of the tracer plume was indeed spatially extended in the transverse and longitudinal direction, we use here for simplicity a pointlike initial particle distribution, localized at the origin at x = 0,

ρ(x) = δ(x). ( 10
)
Regarding the initial velocity distribution p 0 (v), we note that Boggs et al. [12] and Fiori et al. [39] pointed out that more mass entered in the high than in the low conductivity zones, which implies that the initial mass distribution is approximately flux-weighted. In order to set the initial velocity distribution p 0 (v), we assume that the initial plume extension is large enough to be considered ergodic [39]. Under this assumption, p 0 (v) can be set equal to p v (v) given by Eq. ( 2). Note that the model can be conditioned on velocity data within the initial plume if detailed information is available. Below we use this approach to model and predict the longitudinal concentration profiles of the MADE-1 and MADE-2 experiments.

B. Longitudinal concentration profiles

The experimental concentration profiles for the MADE-1 experiment are based on the bromide tracer data from Adams and Gelhar [START_REF] Adams | [END_REF] normalized to the initial mass as given in Fiori et al. [39]. The concentration profiles for the tritium tracer at the MADE-2 experiment are obtained from Salamon et al. [34].

Figures 2 and3 show the experimental data and the model predictions. As explained in Adams and Gelhar [START_REF] Adams | [END_REF] and Salamon et al. [34], the displayed data was obtained by interpolation of concentration values between sampling wells, vertical and transverse integration of the resulting concentration field, and longitudinal averaging over a window of ∆x = 10 m [START_REF] Adams | [END_REF]34] provides the tracer distributions with an averaging window of ∆x = 10 -1 m for MADE-1. We refer to the latter as the fine scale model data.

We first note that the experimental data does not integrate to 1. In fact the mass under the longitudinal profiles for the MADE-1 data integrates to 2.06, 0.99, 0.68, 0.62, 0.54 and 0.43 for the six snapshots at increasing time [START_REF] Adams | [END_REF]39] and to 1.52 and 0.77 at t = 27 d and 328 d for the MADE-2 snapshots [START_REF] Stauffer | Degradation of aromatic hydrocarbons in an aquifer during a field experiment demonstrating the feasibility of remediation by natural attenuation[END_REF]. The explanation of this apparent mass loss has been a matter of debate in the literature. Several authors [30,34,[START_REF] Boggs | [END_REF][START_REF] Molz | [END_REF] attribute the overestimation of mass at early, and underestimation at late times to a measurement bias due to preferential sampling from high conductivity regions and mass transfer into immobile regions. In fact, the analysis of core samples two years after the MADE-1 experiment by Boggs and Adams [START_REF] Boggs | [END_REF] indicate that a substantial amount of mass may have been transferred and stored in immobile zones. Other authors [START_REF] Adams | [END_REF]39] debate that the mass loss could be attributed to a decreasing density of sampling points at distances larger than 20 m downstream from the injection region. This view is supported by the fact that more mass loss is observed for the bromide plume in the MADE-1 [START_REF] Adams | [END_REF] than for the tritium plume in the MADE-2 experiment [START_REF] Stauffer | Degradation of aromatic hydrocarbons in an aquifer during a field experiment demonstrating the feasibility of remediation by natural attenuation[END_REF], for which the density of sampling wells was increased. The tracer distribution in the upscaled TDRW model is by default normalized to 1/θ, which, as outlined above accounts for a possible sampling bias due to adsorption. Figures 2 and3 display model predictions for K g = 5.5×10 -6 m/s, and K g = 10 5 m/s, and K g = 4.3 × 10 -6 m/s at the upper and lower limits of the confidence interval. The model predicts the overall shape of the observed tracer distributions, the slow moving peak, and the development of a forward tail with increasing time. For increasing K g , the velocity spectrum is shifted to higher values. Thus, the peak height and peak retention decrease, the forward tail increases. While the value of K g = 5.5 × 10 -6 m/s provides the best prediction for the peak concentrations, all K g values give qualitatively the same non-Gaussian tracer profiles. This shows the robustness of the salient non-Gaussian transport features predicted by the TDRW approach. The upscaled model accounts for retention of tracer mass in the injection region through persistent low initial flow speeds. The forward tailing is accounted for through the presence of high initial flow velocities in combination with fast tracer motion through particle transitions into spatially persistent fast velocity channels.

The predicted mass in the tails at 126 d to 370 days (MADE-1) is generally higher for x > 25 m than in the experimental data, while reach is relatively well reproduced, see also the corresponding logarithmic plots in Appendix C. As mentioned above, this discrepancy between model and data may be attributed to a decreasing density of sampling points more than 20 m downstream from the injection [START_REF] Adams | [END_REF]39]. The last snapshots at t = 503 d for the MADE-1 and t = 328 d for the MADE-2 experiments show farther reaching forward tails with secondary peaks, which are not captured by the TDRW model. These features in the data may be attributed to a sampling bias.

IV. CONCLUSIONS

We have presented a Lagrangian stochastic model for transport in heterogeneous aquifers that is parameterized in terms of the statistical characteristics of hydraulic conductivity and Eulerian flow speed. The approach is based on a Markov processes for equidistant particle speed, whose steady state distribution is given by the flux-weighted flow speed distribution. The latter is derived from the distribution of the Eulerian flow speed, which is obtained from standardized numerical flow simulations, this means, for unit gradient of hydraulic head and unit geometric mean conductivity. The model is parameterized by the mean, variance and correlation lengths of the logarithm of hydraulic conductivity, advective tortuosity, porosity. These model inputs are transport independent, this means they can be assessed by medium and flow characterization. These features yield a predictive upscaled approach for advective tracer transport in heterogeneous aquifers. Like other stochastic models, this approach relies on the ergodicity of the underlying medium heterogeneity. The model furthermore accounts for linear instantaneous equilibrium sorption.

We use this upscaled model to assess and predict the evolution of the bromide tracer plume of the MADE-1 and the tritium plume of the MADE-2 experiments, which are characterized by strongly non-Gaussian shapes. Relatively low conductivity in the source zone leads to slow peak movement of only a few meters over the duration of the experiments, while part of the tracer moves fast, which gives rise to a pronounced forward tail. The slow moving peak is captured by the upscaled model through the low velocity end in the distribution of initial particle speeds. The developing forward tail is quantified by tracer propagation due to spatially persistent broadly distributed particle speeds. The speed distributions are determined by the medium and Eulerian flow statistics. As mentioned above all model parameters are obtained from transport independent measurements, without recourse to fitting parameters. We find that the purely advective upscaled model provides a robust estimate of the evolution of the overall plume shape. Our results indicate that the large scale transport features of the MADE-1 and MADE-2 plumes may be explained by flow heterogeneity due to spatial variability in hydraulic conductivity. The proposed upscaled model is predictive, transferable to different solutes and hydraulic conditions, and seems to capture the salient heterogeneity mechanisms, and their impact on large scale transport.

Figure 1 .

 1 Figure 1. (Squares) empirical distribution of f (x). The solid orange line denotes the normal distribution with zero mean and variance σ 2 f = 5.9. (Circles) empirical distribution of the logarithm ν(x) of the flow speed. The blue solid line denotes the skewed normal distribution (6).

Figure 2 .

 2 Figure 2. Concentration profiles from (crosses) the MADE-1 data and (circles) the upscaled model for an averaging window of ∆x = 10 m at times (top left to bottom right) t = 49 d, 126 d, 202 d, 279 d, 370 d and 503 d using a pointlike initial distribution. The dashed lines denote the predictions for the values of (green) Kg = 4.3 × 10 -6 m/s and (orange) Kg = 10 5 m/s at the lower and upper limits of the 95 % confidence interval of Kg.

Figure 3 .

 3 Figure 3. Analogous to Figure 2 for the concentration profiles of the MADE-2 experiments at t = 27 d and 328 d.

Figure 4 .

 4 Figure 4. Concentration profiles from (crosses) the MADE data and the upscaled model for an averaging window of (circles) ∆x = 10 m and (solid blue line) ∆x = 10 -1 m at times (top left to bottom right) t = 49 d, 126 d, 202 d, 279 d, 370 d and 503 d.

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Same concentration data as in 2 plotted logarithmically on the concentration axis.

Figure 8 .

 8 Figure 8. (Left panel) Maximum concentrations versus time from the (black crosses) MADE-1 data, and (blue squares) the prediction from the upscaled model with ∆x = 10 m. The dashed lines denote the cumulative distribution of residence times (D1) obtained from the upscaled model for (green) Kg = 4.3 × 10 -6 m/s, (orange) Kg = 10 -5 m/s. The black solid line denotes the prediction for Kg = 6.7 × 10 -6 m/s. (Right panel) Distributions of the corresponding ν = ln(v).

Table I .

 I Geostatistical, and hydraulic data for the MADE site, as well as retardation coefficients for the bromide and tritium tracers of the MADE-1 and MADE-2 experiment according to Bohling et al.[51], Boggs et al.[12], Boggs and Adams

							]
	Value	5.9	5.5 × 10 -6 9.1	1.8	3.6 × 10 -3 0.31	1.2	1

ACKNOWLEDGMENTS

We thank Daniel Fernandez-Garcia for sharing the data for the two concentration snapshots of the MADE-2 experiment, and Alraune Zech and Aldo Fiori for sharing the MADE-1 data with us. The authors acknowledge the financial support of the European Research Council through the project MHetScale (Grant agreement No. 617511), and the Spanish Ministry of Science and Innovation through a Severo Ochoa project (No. CEX2018-000794-S), and the project HydroPore (PID2019-106887GB-C31). Juan Hidalgo acknowledges the support of the Spanish Ministry of Science and Innovation through a Ramón y Cajal fellowship. The experimental data was published originally in Adams and Gelhar [7] and reproduced in Fiori et al. [39], from where the data in the present work is obtained. The TDRW model data is availabe at http://dx.doi.org/10.20350/digitalCSIC/12667.

Appendix A: Flow solution

The hydraulic conductivity fields are generated using the spectral method of Ruan and McLaughlin [59] with an exponential covariance function. The flow problem (3) is solved numerically using a two-point flux finite volume method in 40 realizations of f (x) whose spatial mean and variance are within a 0.05 tolerance interval around the target values of 0 and 5.9. The flow domains have the dimensions 20 h in the direction of the mean gradient, 10 h in the transverse horizontal direction, and 10 v in vertical direction. The resolution is 2 8 × 2 7 × 2 5 voxels.

Appendix B: Numerical model implementation

In the following, we outline the basic steps of the numerical implementation of the stochastic time domain random walk model combined with mass transfer between mobile and immobile regions. First, we note that the steady state Eulerian velocity distribution p e (v) is given by the distribution show in Figure 1. The distribution p v (v) of particle speeds is obtained according to (2) by flux-weighting.

Step 1: Initialization Initial particle velocities v(0) are sampled from p v (v), particle positions x(0) are chosen from the initial particle distribution ρ(x), which here is ρ(x) = δ(x).

Step 2: Propagation of normal scores The initial velocities v(0) are then converted to the initial values w(0) of the normal scores using the map [START_REF] Adams | [END_REF]. The normal scores w(s) are propagated from their initial values w(0) according to the Ornstein-Uhlenbeck process (8), which is discretized by using an Euler scheme. This gives

where ζ(s) denotes a Gaussian random variable with zero mean and unit variance. The discretization ∆s is chosen much smaller than c . Here we set ∆s = 10 -1 m.

Step 3: Propagation of particle position and time The particle positions are incremented at each random walk step by ∆s/χ, the particle times by ∆s/v (s),

The current particle speed v (s) is obtained according to Eq. ( 7).

Steps 2 and 3 are repeated until the maximum simulation time is reached. Particle positions are recorded at steps s if t(s) ≤ t i and t(s + ∆s) > t i with t i the desired observation times.

Appendix C: Longitudinal concentration profiles, Figure 4 shows the longitudinal concentration profiles from the MADE-1 data, and the model prediction in coarse and fine resolution for θ = 1.2. The comparison shows that the averaging window of ∆x = 10 m introduces an artificial broadening of the peak due to oversmoothing, while the tail concentrations are nearly identical. With increasing time, as the peak width increases due to hydrodynamic dispersion, the differences between the coarse and fine scale models decrease.

Figure 5 shows the same data as in Figure 2 plotted logarithmically. The tail contains in general more mass than the experimental data, but reproduces relatively well the overall tendency except for the 503 d snapshot at which mass is recovered at larger downstream distances than predicted by the model.

Appendix D: Predictions for different values of Kg and peak concentrations

Figure 6 compares the model predictions for the MADE-1 data with θ = 1.2 based on the estimated average value of K g = 4.3 × 10 -6 m/s, and the values at the upper and lower limits of the 95% confidence interval. K g = 6.7 × 10 -6 m/s and K g = 10 -5 m/s. For all K g values the model predicts the main plume features of peak localization and forward tail discussed in the main text. The model predictions based on K g = 5.5 × 10 -6 m/s shown in the main text provide slightly better estimates for the concentration peak values at t ≤ 370 d than the ones using the mean K g = 6.7 × 10 -6 m/s.

The concentration peak is localized at all snapshots at x = 5 m, which indicates that it does not move beyond the averaging window between x = 0 m and 10 m. This is supported by the fine scale prediction. The upscaled model predicts peak heights and widths that are qualitatively and quantitatively similar to the experimental observations. The peak concentrations from the experimental data, and the corresponding model predictions are displayed in Figure 8. The upscaled model predicts the overall decay of the peak concentration. The decay behavior in the localized peak is in fact a measure for the residence time in the injection region. The mass m 0 (t) remaining in the region between 0 m and 10 m is equal to the probability that the residence time τ 0 of a particle is larger than t,

where ψ 0 (t) denotes the residence time distribution. As long as the mass in the injection region of size ∆x is larger than the mass in all other bins, the peak concentration is equal to c 0 (t) = m 0 (t)/∆x. Figure 8 shows the evolution of c 0 (t) from experimental data and model predictions obtained from the upscaled model. In order to highlight the impact of the value of geometric conductivity, we also plot the decay of the maximum concentration for the values K g = 6.7 × 10 -6 m/s and the values at the upper and lower limits of the 95% confidence interval, K g = 10 -5 m/s and K g = 4.3 × 10 -6 m/s, as well as the corresponding distributions of ν = ln(v). The former underestimates the observed maximum concentrations because, as shown in the right panel of Figure 8 the distribution of ν is characterized by higher values than for the lower K g . This feature together with the persistence represented by the correlation length c of particle velocities leads to a relatively fast tracer release from the injection domain compared to the other values for K g . The value of K g = 4.3 × 10 -6 m/s overestimates the maximum concentration because of strong tracer retention at the injection region due to low initial speeds. Nevertheless, the predictions using these K g values provide robust prediction of the salient non-Gaussian behaviors of the tracer profiles.


Hippolyte D'albis

Ikpidi Badji
email: ikpidibadji@gmail.com

Najat El Mekkaoui
email: najat.el-mekkaoui@dauphine.fr

Julien Navaux
email: jnavaux@uottawa.ca

Private asset income in France: Is there a breakdown of intergenerational equity between 1979 and 2011?

Keywords: National Transfer Accounts, Private Asset Income, Intergenerational Equity. JEL codes. E01, E21, D14

We use the National Transfer Accounts methodology to calculate private asset income by age for the years 1979-2011. We analyze age profiles using three indicators of intergenerational equity.

Monetary asset income shows no evidence of generational breaks to the benefit of the baby-boom generation. On the contrary, baby-boomers suffered from the high interest rates that they paid to become homeowners. Imputed rents show an obvious breakdown of intergenerational equity when we use an inter-age and intergenerational indicator. This indicator compares the per capita asset income at a given age with the average asset income of people aged 18-85. It gives the relative situation of one age group compared to its contemporaries and it also gives the relative situation of one generation when we compare birth cohorts over time. We find that the cohort born in 1950 benefited from a better position than their successors. Moreover, the cohorts born before the war and during the war appear to be even more favored than the baby-boomers. The cohorts born in 1930 and in 1940 have a better situation than the previous generations and a better position than the following generations.

Introduction

The French PAYG system reached a peak in 2014. That year, public pensions accounted for 14.1% of GDP (Conseil d'Orientation des Retraites 5 , 2017). From 2014 forward, that share should progressively fall, provided that the GDP growth rate maintains an acceptable level. If the GDP growth rate remains stable at 1.3%, pensions could account for 13.3% of GDP by 2050. A higher growth rate would result in a smaller share of GDP -for example, with a stable 1.8% growth rate, pensions as a share of GDP could fall to 12.1%. Contrary to popular belief, France has in fact implemented pension reforms to address the issue of population aging. Although the reforms were justified on grounds of sustainability, they will be costly for the pensioners. Replacement rates will decrease strongly from the 1950's birth cohorts.

Whereas a 68-year-old born in 1947 could count on a replacement rate of 54.3%, a 68-year-old born in 1980 will see his replacement rate fall to 47.5%6 . Consequently, future pensioners will have to find a solution to the conundrum: how to maintain its consumption level over time as public pensions decline.

Should they rely on their families or on themselves? Family solidarity with transfers from younger generations to older generations is not a realistic solution in France. Nowadays, the direction of family transfers is the opposite. French pensioners transfer a net amount equivalent to 0.9% of GDP to younger generations7 . However, given that pensioners concentrate more assets than younger generations, another solution may come from retirees' assets.

In 2011, private asset income accounted for 17.1% of GDP, i.e. 3.6 points higher than public pensions for the same year, and the share of private asset income received by persons aged 60+ was 47.5%. However, French pensioners suffered from the low interest rates implied by the monetary policy of the European Central Bank. In 2016, long-term interest rates were lower than 0.47 % in France, while they were higher than 1.8% in the USA.

In a nutshell, future pensioners could suffer from pension reforms without any guarantee of benefiting from higher rates of return on their assets. It should be pointed out that the consumption level of pensioners is not restricted to monetary flow. It also includes imputed rents that are a significant part of private asset income. In 2011, rents accounted for 30.3% of total private asset income. Moreover, the share of homeownership increased significantly for people aged 50 and over 3 between 1986 and 2010 [START_REF] Arrondel | Les comportements patrimoniaux des ménages en France: Évolutions et déterminants entre 2004 et 2010[END_REF]. This increase was even greater for people aged 70 and over. People who bought their first home before the price boom that started in the late 1970s also benefited from a strong leverage effect.

We expect that the combined effects of pension reforms, housing prices, and the economic climate will strongly modify the allocation of asset income across generations. In France, the most recent studies on the allocation of consumption and labor income (d'Albis et al., 2017;[START_REF] D'albis | Intergenerational inequalities in standards of living in France[END_REF], on public transfers [START_REF] Chojnicki | L'impact des réformes des retraites sur l'équité entre générations[END_REF][START_REF] Navaux | Intergenerational transfers in France: Stable trends and breaks between age groups[END_REF], and on financial wealth [START_REF] Bernard | French households financial wealth: which changes in 20 years?[END_REF] showed no evidence of an intergenerational break to the benefit of the baby-boomers. However, these studies do not address the issue of the allocation of private asset income between age groups and generations.

Our study is focused on all these aspects of asset income. It deals with economic flows and with imputed rents. We use the National Transfer Account (NTA) methodology to investigate the issue of intergenerational equity. We introduce age profiles for financial and real estate income, household interest rates, capital income of self-employed individuals, and for imputed rents. The French NTA project has been conducted with the aim of measuring time series over the longest period, i.e. between 1979 and 2011. On this note, France has now joined the small group of countries that have constructed NTA time series, such as the United States [START_REF] Donehower | The changing shape of the economic lifecycle in the United States, 1960 to 2003[END_REF], Sweden [START_REF] Lindh | Back-casting National Transfer Accounts in Sweden from 1800 to 2009[END_REF], Taiwan [START_REF] Lai | Who supports the elderly? The changing economic lifecycle reallocation in Taiwan, 1985 and 2005[END_REF], and Australia [START_REF] Temple | Mature Age Labour Force Participation and the Life Cycle Deficit in Australia: 1981-82 to 2009-10[END_REF]. Access to long term time series for French age profiles is a necessary prerequisite for comparing several generations at the same age.

We introduce three indicators of intergenerational equity with the aim of obtaining a broad analysis of intergenerational differences. The first indicator describes the relative situation of an age group with respect to the situation of another age group over time. Our second indicator compares the standard of living of successive generations at the same age. The third indicator is a combination of the two previous criteria. It is both an inter-age indicator and also an intergenerational indicator. It computes relative measures of income by dividing net resources at age and date by their cross-sectional average. This paper is organized as follows. The first section describes NTA methodology and the data sources for the age profiles. The second section introduces the three intergenerational equity indicators that we selected to analyze the NTA profiles. The third and the fourth sections present the empirical results for monetary asset income (section 3) and for imputed rents (section 4) and the last section concludes.

Methodology and data

NTA is an international project whose purpose is to "improve understanding of how population growth and changing population age structure influence economic growth, gender and generational equity, public finances, and other important features of the macro-economy" notably by creating National Accounts broken down by age. It computes age profiles using survey and administrative data that is readjusted by using aggregate controls drawn from National Accounts. The theoretical basis of the methodology is available in [START_REF] Lee | Population Age Structure, Intergenerational Transfer, and Wealth: A New Approach, with Applications to the United States[END_REF]Mason (2011a, 2011b) and the detailed methodology for building those accounts is explained in a manual edited by the United Nations in 2013. Some recent developments using NTA are presented in d' Albis andMoosa (2015), d'Albis et al. (2017) and [START_REF] Temple | Mature Age Labour Force Participation and the Life Cycle Deficit in Australia: 1981-82 to 2009-10[END_REF].

Private asset income is broadly broken down between property income and capital income. Table 1 presents the detailed breakdown of the aggregate and gives the amount of the flows for 2011. In this article, we focus on income received by households and exclude those received by corporations and nonprofit institutions serving households (NPISHs). The computation of our two main variables, based on households' private property income and private capital income, is obtained as follows. Private property income of households includes financial asset income, which amounted to 113.2 billion euros in 2011, and real estate income (including those obtained from land), which amounted to 39.7 billion euros.

Moreover, interest payments on loans are subtracted. In 2011, households gave 25 billion euros to corporations for that purpose. Private capital income of households includes two groups of income inflows: the private capital share of mixed income8 , which represented 27.5 billion euros, and private capital income on owner-occupied housing, also referred to as imputed rents, which amounted to 106.6 billion euros. Thus, private capital incomes include both real flows (financial asset income, real estate income, interest, capital share of mixed income, and fictive flows (imputed rents). In what follows, we aggregate real flows in a single variable that we have called "monetary asset income", corresponding to an aggregate net amount of 155.4 billion euros in 2011, and we analyze the 106.6 billion euros of imputed rents separately.

Two different surveys are used to calculate the age profiles: the French Household Expenditures survey (Budget de famille) and The French Wealth survey (Actifs financiers and Patrimoine). Financial asset income, interests, the capital share of mixed income and imputed rents are taken from the Household Expenditure surveys conducted in 1978-1979, 1984-1985, 1989, 1994-1995, 2000-2001, 2005-2006 and 2010-2011. Real estate incomes are taken from the Wealth surveys conducted in 1986, 1991, 1997-1998, 2003-2004 and 2009-2010. It should be mentioned that private income is reported at the household level (except mixed income that is available for each individual) 9 and the surveys give no information about the owner of each asset.

According to NTA methodology, asset income inflows should be assigned to the household head. In order to provide a better individualization within the household, we choose to share each asset income equally among spouses. Unfortunately, the 1979 to 1995 expenditure surveys do not provide figures for imputed rent. We had to estimate them using the characteristics of housing. The procedure is similar to the one used in (d'Albis et al., 2015(d'Albis et al., , 2017)). The imputed rent of homeowners is calculated using the following equation:

Rent i = exp (X i β + residual i )
where X i is the vector of the variables (region, urban units, surface area, number of rooms, etc.) of the rent equation for observation i and where β is the vector of the estimated coefficients of the rent equation. This regression is conducted separately for houses and for apartments. In order to obtain correct rent distribution, the imputed residual must have the same distribution as the residuals taken from the rent equation. As the rent equation residuals are heteroscedastic and non-Gaussian, they cannot be expressed as a normal distribution. The appropriate residual imputation method is the Hot Deck method, which involves randomly selecting an estimated residual using the estimation from the rent equation. This residual is then imputed to housing "similar" to the one from which we selected the estimation residual and for which we have to calculate the imputed rent.

Once the average age profiles are calculated, we smooth each profile by using the [START_REF] Friedman | A variable span scatterplot smoother[END_REF] method recommended in the NTA methodology (Appendix B of United Nations, 2013). We highlight the fact that each survey provides distant waves from several years and that Expenditure and Wealth surveys were not conducted during the same years. Consequently, we chose to interpolate the variables for all ages using polynomial functions. One should note that the wealth surveys are not available for 1979 and 1985 or for 2011. Thus, it is impossible to use the polynomial smoother for the first and the last year of real estate income because the interpolation methodology requires obtaining observations between the years that we interpolate. 10 We assume that the shape of the age profiles between 1979 and 1985 are equivalent to the shape of the year 1986 and the shape of the profile for the year 2011 is equivalent to that of 2010. Finally, we obtain NTA age profiles for each variable over a period of 32 years, from 1979 to 2011.

Figure 1

Three indicators of intergenerational equity

We selected three indicators of intergenerational equity with the aim of obtaining a broad analysis of intergenerational differences.

The first indicator describes the relative situation of retirees with respect to the situation of younger people. The evolution of this indicator over time is a way to evaluate intergenerational equity: in particular, a constant ratio can be interpreted as an indication of equity between age groups. This indicator provides a useful instantaneous picture, which describes the relative situation of two age groups at a given time [START_REF] Blanchet | Mesurer l'équité intergénérationnelle: Le choix des indicateurs[END_REF]. For instance, [START_REF] Navaux | Intergenerational transfers in France: Stable trends and breaks between age groups[END_REF] showed that in 2011, per capita public transfer inflows of French people aged 60 to 85 were 3.1 times higher than per capita public transfer inflows of French people aged 25 to 59. The analysis of fluctuations of this ratio over time is also interesting. The baseline situation may be based on a stable age-to-age ratio as observed by [START_REF] Becker | The Family and the State[END_REF] The second indicator compares the asset income of successive generations at the same age. For this indicator, the standard of intergenerational equity requires that, at each age, every generation's asset income has to be at least equivalent to that of the previous generation [START_REF] Bonnet | Un inévitable conflit des générations ?[END_REF]. For instance, we ask ourselves whether the situation of the baby-boomers is better than the situation of the following generations at the same age. If this were the case, it would represent a break in intergenerational equity to the benefit of baby-boomers. However, this indicator is a minimum standard for intergenerational equity, because a period of economic growth should result mechanically in an improvement in the standard of living from one generation to another, as it had been shown in France for labor income, consumption, and public transfers [START_REF] Navaux | Intergenerational transfers in France: Stable trends and breaks between age groups[END_REF]. This indicator must be completed by a more ambitious one.

The third indicator is a combination of the two previous criteria. It is, at the same time an inter-age indicator and an intergenerational indicator. It computes relative measures of asset income by dividing net resources at age and date by their cross-sectional average [START_REF] Blanchet | Pension and Generational Histories in a Simple Demo-economic Model: the Case of France[END_REF]. In this paper, we divide per capita asset income of persons aged at time with the average asset income of people aged 18 to 85:

Per capita current asset income of persons aged at time Average of per capita asset income by age unweighted by demography of persons aged 18 85 at time By calculating this ratio for each year, we are able to follow each generation over the 32-year time window . This is an inter-age indicator, because for a given year, the indicator reveals what age group is in the best position in comparison with the other age groups. This is also an intergenerational indicator because it is possible to determine whether a given birth cohort maintains a better position than its contemporaries over time. A golden generation would be a birth cohort that holds the best position over its entire lifecycle. This indicator is stronger than the second indicator because it is growth neutral in the sense that a generation can maintain a favored position over time even with high economic growth. It was not the case for the second indicator that depends mainly on GDP growth.

The three indicators that we have selected do not cover all the possible ways to study intergenerational equity but provide a good coverage of the issue. Moreover, the period of 32 years that is available for France is not enough to construct age profiles for the complete life cycle of generations. Consequently, we cannot measure intergenerational equity indicators based on a balance sheet of several generations over their life cycle, as it would require overlapping generation models [START_REF] Chojnicki | L'impact des réformes des retraites sur l'équité entre générations[END_REF] or generational accounting exercises [START_REF] Kotlikoff | Generational Accounting. Knowing Who Pays, and When, for What We Spend[END_REF]. One way to overcome the difficulty is to build pseudo-cohorts as in d'Albis and Badji (2017). Figure 2(a) shows that the per capita normalized age profile remained stable over time for people between 40 and 75 years of age. However, the age profile diverged for younger and older people. After 75, the normalized age profile increased between 1979 and 1984, decreased between 1984 and 2000, and remained stable afterward. A symmetric evolution was observed for younger people between 18 and 40 years of age. Monetary asset income became strongly negative between 1979 and 1984, then increased up to 1995, and remained quite stable thereafter. These fluctuations can be partly explained by the interest paid by individuals, most notably before 40 years of age. For instance, real interest rates on mortgages increased considerably between 1979 and 1985, then remained stable until 1992, when they began to decrease [START_REF] Wilhelm | L'évolution actuelle du crédit à l'habitat en France est-elle soutenable ?[END_REF].

No large differences in monetary asset income between generations

The impact of the high interest rates in the 1980s also appears in Figure 2(b) which shows the ratio of per capita monetary asset income of people aged 60 to 85 to per capita monetary asset income of people aged 18 to 59. The ratio reaches a peak in 1984 because of the high interest rates paid by young people.

In 1984, people aged 60 to 85 received 7.2 times the average monetary asset income of people aged 18 to 59. Thereafter, the ratio follows an L-shaped curve with the lowest point of 2.3 in 2003. A second peak, of 3.2, is reached in 2008 and the ratio decreases slightly again after 2008. The high peak of 1984 draws our attention as it could indicate an imbalance between generations. In particular, one could ask whether a generation could have a lower monetary asset income than previous generations at the same age. Figure 2(c) describes the evolution of monetary asset income by birth cohort in constant euros. We assume that the evolution of monetary asset income is equitable for generations since a generation has a higher income than the previous generations at the same age. On the whole, Figure 2(c) shows that monetary asset income increases from one generation to another between 1979 and 2011. For instance, at age 40, people who were born in 1960 get more than two times the average constant monetary asset income of people who were born in 1950. At age 60, monetary asset income of people who were born in 1940 was 47% higher than people who were born in 1930 at the same age. D' Albis et al. (2017) showed that successive generations benefited from economic growth through labor income, and our results

show that it was also the case for monetary asset income. Moreover, the figure shows a specificity for the 1950 and 1960 birth cohorts because these are the only generations who received a negative monetary asset income at young ages. The 1960 birth cohort reaches a minimum of -203 constant euros in 1985 at age 25. This effect is even stronger for the 1950 birth cohort who receives a negative monetary asset income of -716 constant euros in 1982 at age 32. We evidence once again the impact of high interest rates in the mid-1980s that did impact these generations just when they bought their first home. Figure 2(c) shows also that people who were born in 1980 had systematically lower asset income than people who were born in 1970 when we compare them at the same age, i.e. up to age 31. The window of observation is rather small but it shows that this generation experienced a deteriorated situation at the beginning of its working life.

Significant generational effects for imputed rents

The normalized age profiles presented in Figure 3(a) show that imputed rents follow two successive evolutions.12 First, imputed rents decrease over age after age 50 in 1979, in 1984, in 1989, and in 1995. For these four survey waves, the slope is less and less steep. In 2000, in 2005, and in 2011, the age profile increases after age 50. This shift at older age could be explained by a significant increase in the ownership rate for people aged 60 and more after This increase is stronger after age 70 than between 60 and 70 [START_REF] Bonvalet | Les effets de la politique du logement sur l'évolution du taux de propriétaires en France[END_REF]. We find a similar dynamic in Figure 3(b), which represents the ratio of per capita imputed rents of people aged 60 to 85 to per capita imputed rents of people aged 18 to 59. In 1979, the per capita imputed rents of people aged 60 and over was equal to the per capita imputed rents of younger people. From 1979 to 2000, the ratio increased to 1.5 in 2000, which means that the per capita imputed rents of people aged 60 and over was 50% greater than the per capita imputed rents of people aged 18 to 59. Thereafter, the ratio decreases again to 1.3 in 2011. The two successive movements of figure 2(b) can be explained by the situation of older individuals on the housing market, as we showed in figure 3(a), but also by the macroeconomic environment that impact the young. From the late 1970s to the late 1990s, young generations faced a substantial increase in housing prices [START_REF] Friggit | Quelles perspectives pour le prix des logements après son envolée?[END_REF], high unemployment rates, and interest rates that were still high [START_REF] Clerc | Les inégalités entre générations depuis le baby-boom[END_REF]. Consequently, young people endured a difficult macroeconomic environment that impacted their capacity to become homeowners [START_REF] D'albis | Evolutions démographiques et marché de l'immobilier neuf[END_REF]. After 2000, the share of homeowners rose sharply for young people (especially between 25 and 35 years of age) because of low interest rates [START_REF] Clerc | Les inégalités entre générations depuis le baby-boom[END_REF], better income expectations [START_REF] Gobillon | L'impact des contraintes d'emprunt sur la mobilité résidentielle et les choix entre location et propriété[END_REF], and positive profit expectations due to increasing housing prices [START_REF] Poterba | Taxation and Housing: Old Questions, New Answers[END_REF]. The share of homeowners also increased for people aged 60 and over but was not sufficient to compensate for the effect of interest rates on young people.

Fluctuations in the macroeconomic environment have an impact on imputed rents of generations and thus call into question intergenerational equity. Figure 3 Moreover, for most of the years, those generations have a better situation than the previous generations and the following generations when we compare them at the same age. It appears that these generations occupied an even better situation than baby-boomers born in 1950. In fact, people born in 1950 have a better situation than the cohorts born in 1960 and 1970 at each age15 . The difference seems to be significant between the baby-boomers born in 1950 and younger people at the same age, but it also significant between the 1950 cohort and the previous generations. The generations born between 1960 and 1980 have a similar situation at the same age, which demonstrates a relative equity between them. Thus, the cohort born in 1950 indeed benefited from a favored position compared to their successors but also had a worse situation between ages 35 to 50 compared to cohorts born in 1940 and 1930, which appear to be the most favored generations. This relative situation on the housing market does not imply a higher standard of living as measured by consumption or revenues including imputed rents (d'Albis and Badji, 2017). generations that followed. This is due to the prevailing impact of interest rates that exceeds the impact of imputed rents for the 1950 cohort. On the whole, this cohort received a lower asset income than the 1960 cohort, the 1970 cohort, and the 1980 cohort at each age between 1979 and 2011. Figure 4 also reveals that monetary asset income tends to hide all the other generational breakdowns that were evidenced for imputed rents when we consider total asset income.

Conclusion

In France, several studies have been conducted recently on intergenerational equity. They demonstrate that baby-boomers were not favored with regard to the allocation of consumption, of labor income, of public transfers, and of financial assets (d'Albis et al., 2017;[START_REF] D'albis | Intergenerational inequalities in standards of living in France[END_REF][START_REF] Chojnicki | L'impact des réformes des retraites sur l'équité entre générations[END_REF][START_REF] Navaux | Intergenerational transfers in France: Stable trends and breaks between age groups[END_REF][START_REF] Bernard | French households financial wealth: which changes in 20 years?[END_REF]. Our paper extends this conclusion to monetary asset income.

We show that baby-boomers didn't benefit from a favored position in comparison with the other generations, because they suffered from the high interest rates of the 1980's. The cohort born in 1950 received strong negative monetary asset income when they were young adults, while the next generations had slightly negative or slightly positive monetary asset income at the same age.

On the contrary, imputed rents imply several breaks of intergenerational equity when we follow the relative position of each generation with respect to their contemporaries over time. First, the cohort born high unemployment rates and high interest rates. However, it does not mean that the cohort born in 1950 is a golden generation, because the generations born in 1930 and in 1940 appear to be much more favored than the baby-boomers. They benefited from a better position than the previous generations and the following generations when we compare them at the same age.

Although there is no evidence of any generational break to the benefit of baby-boomers for actual economic flows, our results for imputed rents reveal an important issue about wealth. We know from [START_REF] Piketty | Capital in the 21st century[END_REF] that the share of wealth in National Income and the share of wealth transfer in National Income increased significantly since World War II and that it might strengthen over the coming years. Therefore, our current research focuses on an extension of National Transfer Accounts to wealth accounts, following the primary works of [START_REF] Lee | Population Age Structure, Intergenerational Transfer, and Wealth: A New Approach, with Applications to the United States[END_REF]. We will also carry out an update of economic flow accounts up to 2017 with the next Budget Survey that should be available in 2018. Moreover, asset income is highly concentrated on top income deciles. This process is reinforced since wealth is more and more concentrated, particularly because of the growing economic weight of inheritance [START_REF] Piketty | Capital in the 21st century[END_REF].

Consequently, National Transfer Accounts by income level are necessary to assess whether birth cohorts of each income decile improved their situation with respect to the previous generations that were in the same income decile. It is nevertheless necessary to create a new NTA methodology since macroaggregates are not available for income deciles.

Figure 2

 2 Figure 2 introduces the per capita age profiles normalized to 1 in Figure 2(a), which are computed such that the surface below the curve is equal to 1, as well as the three intergenerational equity indicators in Figure 2(b), Figure 2(c), and Figure 2(d).

  Figure 2(d) is useful to assess if a generation has a better position than its contemporaries and how its position evolves over time. Three generational effects are revealed between 1979 and 2011. First, there exists a breakdown of intergenerational equity for the cohorts 1910 to 1930. In fact, the 1910 cohort had a worse position than the 1900 cohort at each age, and the 1930 and 1920 cohorts had a worse position than the 1910 cohort at each age. 11 Secondly, between age 40 years and age 70, all the cohorts have the same relative position. It does evidence a fair distribution of monetary asset income for these age groups. Thirdly, the negative impact of high interest rates of the mid-1980s evidenced by figure 2 (a), figure 2(b) and figure 2(c) is also very strong in figure 2(d). The 1950 birth cohort and the 1960 birth cohort had a worse position than the 1970 and the 1980 cohorts when they were young. After 40, the two cohorts recovered a similar position to previous generations when we compare them at the same age.In a nutshell, monetary asset income involves generational breaks between 1979 and 2011 but the babyboom generations did not benefit from their position over this period. On the contrary, the 1950 cohort experienced a difficult position when they were young by comparison with the 1960 and the 1970 cohorts. Finally, one should keep in mind the limitations of real estate income age profiles, since no 12 surveys are available between 1979 and 1985 and for 2011, and because real estate income represents a significant part of monetary asset income (25.5% in 2011, see table1).

Figure

  Figure 2: Per capita monetary asset income between 1979 and 2011

  (c), which represents cohort profiles normalized in constant euros 13 , shows that imputed rents increase from one generation to another at the same age, except from 2009 to 2011, when people aged 48 and over had lower imputed rents from one generation to another. First, imputed rents increased sharply between 1979 and 2009 for people born between 1900 and 1950. When we compare them at the same age, imputed rents increased by 223% between the 1900 and 1910 cohorts, they increased by 165% between the 1920 and 1930 cohorts, and by 109% between the 1940 and the 1950 cohorts 14 . Imputed rents increased slowly from the 1950 cohort to the 1960 cohort, and even more slowly between the 1960 cohort and the 1970 cohort. Imputed rents stagnated between the 1970 cohort and the 1980 cohort: they increased only by 2.97% in constant euros between 18 and 31 years of age.

Figure 3

 3 Figure 3(d) indicates whether a generation systematically holds a favored position for the entire period compared to other generations at the same ages. We identified several generational breaks between 1979 and 2011. According to our third intergenerational equity indicator, the 1910 and 1930 cohorts benefit from a better situation than the generation born ten years before, except for certain ages. Thus, there is a clear generational break at the expense of older generations. The cohorts born in 1930 and 1940 seem to have a similar situation because their curves almost overlap between ages 49 and 71.

Figure

  Figure 3: Per capita imputed rent income between 1979 and 2011

Figure

  Figure 4: Per capita private asset income between 1979 and 2011 Ind. 3: Inter-age and intergenerational indicator*

  holds a better position than the cohorts born in 1960 and in 1970. It reflects the high housing prices that young generations faced from the late 1970's to the late 1990's, combined with

Table 1 : Asset income in the NTA methodology -Aggregates for France, 2011 (billion €)

 1 Private property income received by corporations and NPISHs are generated from assets owned by corporations and NPISHs while private capital income of corporations and NPISHs are generated by corporations and NPISHs themselves.

	Private asset income	=	Private property income	+	Private capital income
	352.1		63		289.1
			=		=
		Private property income inflows	-	Private property income outflows	Private capital income on owner-occupied housing
		667.4		604.4	106.6
		=		=	+
		Financial asset income		Household interest expense	Private capital income from mixed income
		113.2		25	27.5
		+		+	+
		Real estate income		Other private property income outflows *	Private capital income of corporations and NPISHs **
		39.7		579.4	155
		+
		Private property
		income received
		by corporations
		and NPISHs **
		514.5
	Source: authors' calculations, National Accounts (INSEE). Flows in red are those that are not considered in this article.
	* Other private property income outflows include interests paid by corporations and NPISHs, dividend outflows from the
	corporations, royalties, and rents for land (out of housing rents) and deposits paid by non-financial corporations and
	households. Thus, a very small share of household private asset income is included in "other private property income outflows"
	but it cannot be attributed to individuals.
	**

Breakdown of per capita private asset income in 2011

represents components of private asset income by age in 2011. Financial asset income rises up to 4,771 euros per capita and per year at 69 years of age and decreases slightly to 3,635 euros at 85 years of age. It remains stable thereafter. Real asset income increases up to age 66, the age at which it reaches a peak of 1,471 euros, and remains stable after age 80. The stability of financial asset income after 85 and of real asset income after 80 is not explained by the real shape of the age profiles. In fact, given the small number of observations at old age, the quality of the age profiles extracted from the surveys is very poor. Consequently, for each year we assume that the value of financial asset income after age 85 is equal to the value at age 85 and the value of real asset income after age 80 is equal to the value at age 80. The age profiles of the capital share of mixed income and of interest do not exhibit the same limitations because they are concentrated on working age groups, for which we can obtain many observations. More precisely, the shape of mixed income follows an "M" curve. The age profile reaches a first peak of 858 euros at age 47 and reaches a second higher peak of 1,246 euros at age 60. Imputed rents rise from the beginning of working life up to 43 years of age when it reaches a peak of 2,651 euros per year. Thereafter, imputed rents remain quite stable, as they fluctuate between 2,491 and 2,792 euros. Finally, interest income is greater than 500 euros between ages 29 and 55. Figure 1:

  for public transfer inflows between 1940 and 1983 in the United States and as observed by d'Albis and Navaux (2016) between 1979 and 2011 in France. In this article, we compare per capita by retired people in France (aged 60 and over) that we would like to compare with adults (aged 18 to 59). Due to the bad quality of data after age 85, we defined our indicator as follows:

	8 10 asset income received 60 85 18 59 ,		8 10
	2 4 6 We therefore consider only observations before age 85. Thousand euros				2 4 6
	0									0
	0	10	20	30	40	50	60	70	80	90+
	-2									-2
					Age
		Real estate income		Financial asset income
		Capital share of mixed income	Imputed rents
		Interest income

Institution that advises the French prime minister for the public pension system.

With annual productivity gains of 1.3%.

Data from the National Transfer Accounts for France in 2011.

In the NTA methodology, the mixed income received by self-employed individuals is allocated between labor and capital. The NTA manual (United Nations, 2013) suggests separating the share of capital and labor in total earnings with a 1/3-2/3 rule in order to calculate the aggregate of each component. In France, we made a more precise evaluation of those aggregates using the method proposed by[START_REF] Askenazy | Le partage de la valeur ajoutée[END_REF]. Our detailed methodology is presented in d'Albis et al. (2015).

The shape of the age profile of the capital share of mixed income is equivalent to the shape of the labor share of mixed income, see d'Albis et al. (2015'Albis et al. ( , 2017)).

The same method was applied in France for the life cycle deficit(d'Albis et al., 2015(d'Albis et al., , 2017)).

National Transfer Accounts include per capita imputed rents by age. Unfortunately, it does not provide any information about the price versus quality effects, which can both impact rent growth[START_REF] Friggit | Mission sur les aides personnelles au logement[END_REF].

It would have been better to normalize the cohort profiles of Figure3(c) by using the index of the rents actually paid by tenants. INSEE produces this index, but only since 1990. Consequently, we chose to normalize the cohort profiles in constant euros with the Consumer Price Index (approximately 6% of this index is composed of the rents actually paid).

We compare cohorts in the same age groups, depending on the available data between 1979 and 2011. We compare the cohorts 1900 and 1910 after 79 years old, the cohorts 1920 and 1930 between 59 and 81 years old and the cohorts 1940 and 1950 between 39 and 61 years old.

Except between the cohorts 1950 and 1960 at 49 years old and 50 years old.


Keywords: Electromagnetic analysis, magnetic flux, magnetic forces, optimization, reluctance network (RN)

The aim of this paper is to propose an efficient and unique reluctance network (RN) model for all translator displacements of a linear hybrid motor. This model has to consider several translator positions to analyze the linear motor performances. The proposed RN model takes into account the magnetic characteristic and the flux leakage. Consequently, the proposed model gives precise results of the electromagnetic characteristic. The RN model considers an interesting solution to deal quickly with an optimization task of a large number of constraints and parameters. The developed model is coupled with sequential quadrating programming method for an optimization process of the linear machine. The obtained results confirm the efficiency of this method for modeling and optimization of the linear machine.

T known as flux tubes. Each flux tube is associated with a magnetic reluctance, which depends on circuit dimension and the magnetic behavior. The originality of this model consists HE LINEAR hybrid motors (LHMs) consist of an efficient solution for precise linear displacements applications [START_REF] Krishnan | Switched Reluctance Motor Drives: Modeling, Simulation, Analysis, Design, and Applications[END_REF]. Many research papers work on the design and the optimization of these machines [START_REF] Szabó | A novel double salient permanent magnet linear motor[END_REF]. The finite element method (FEM) is generally used when precise results are indispensable [START_REF] Du Peloux | Automatic generation of sizing static models based on reluctance networks for the optimization of electromagnetic devices[END_REF]. However, this method needs large computation time, limiting the number of parameters that can be taken into account. Moreover, the FEM cannot be coupled with rapid optimization algorithms [START_REF] Raminosoa | Time-stepping simulation of synchronous reluctance motors using a nonlinear reluctance network method[END_REF]. The second method used in the design of linear machines is based on the reluctance network (RN) model [START_REF] Du Peloux | Automatic generation of sizing static models based on reluctance networks for the optimization of electromagnetic devices[END_REF]- [START_REF] Nguyen-Xuan | Efficient reluctance network formulation for electrical machine design using optimization[END_REF]. In this paper, we are concerned with the development of a new RN model for modeling and optimization of linear machines. A 3-D view of the studied machine is shown in Fig. 1. The main dimensions of the proposed LHM are shown in Fig. 2.

The double-sided LHM mechanical dimensions are summarized in Table I.

In Section II, a unique RN model for all translator positions is presented and discussed. The comparison of results with that of the FEM is presented in Section III. The optimization of the machine is presented and discussed in Section IV. The optimization task consists in reducing the global mass of the machine. The electromagnetic characteristics of the optimized structure are also presented and verified using the FEM. The conclusions and perspectives are given in Section V.

II. UNIQUE RELUCTANCE NETWORK MODEL

FOR ALL DISPLACEMENTS The concept is based on the subdivision of the magnetic circuit of the machine into a sufficient number of elements in its ability to calculate the electromagnetic characteristics for various positions from unaligned to aligned positions using a simple run and a unique topology of model. An elementary structure of the LHM with the respective reluctances is shown in Fig. 3.

A. Reluctances of the Air Gap

The air gap tubes are composed of straight lines with a constant length. However, the section of the air gap increases with the displacement x. Hence, the air gap reluctance is parameterized with the displacement, and it is given by the following expression:

e = e μ 0 • x • Tt . ( 1
)
The salient in the mobile canalizes a non-neglected part of the magnetic flux. Hence, the salient is modeled by reluctances expressed as follows:

s = Tpl + e μ 0 • (Tw -x) • Tt . ( 2
)

B. Reluctances of the Stator

Generally, the magnetic reluctance of a ferromagnetic material, with a length l and a section S, can be written as

fer = l μ 0 • μ r • S ( 3
)
μ r is the relative permeability, which depends on the saturation of the ferromagnetic material. Hence, the ferromagnetic circuit of stator is divided into three different magnetic reluctances connected in series. These reluctances depend on both the stator geometrical parameters and the magnetic permeability.

C. Reluctances of the Permanent Magnet

Each permanent magnet inserted in the stator is modeled by a magnetomotive force and a reluctance parametrized with the geometrical and magnetic parameters of the permanent magnet. This magnetomotive force and the reluctance of the permanent magnet can be written, respectively, as

E a = B r μ 0 • μ a L a (4) a = 1 μ 0 • μ a L a S a (5)
where μ a is the relative magnetic permeability of the permanent magnet. L a and S a are, respectively, the length and the section.

D. Reluctances of the Translator Part

The translator part is moved from unaligned to aligned positions. Consequently, the translator canalizes different flux tubes, which depend on the translator position. The RN model developed must be available for each translator position. Hence, the translator part is modeled by different reluctances. The flux in the second translator teeth, which must be aligned with those of the stator in the aligned position, increases with the displacement. Hence the reluctance of the teeth is inversely proportional to the displacement, and it can be written as

21 = Tpl μ 0 • μ r • x • Tt . ( 6
)
Applying the same methodology, the reluctances of other parts are determined. The reluctances 11 and 31 are, respectively, given as

11 = Tpl μ 0 • μ r • Tw • Tt (7) 31 = Tpl μ 0 • μ r • (Tw-x) 2 • Tt . ( 8
)
The reluctances of the cylinder head depend on both the teeth position and the displacement. The reluctances 12 , 13 and 22 are, respectively, given as

12 = Tch μ 0 • μ r • Tw • Tt (9) 13 = Tch μ 0 • μ r • (Tw-x) 2 • Tt (10) 22 = Tch μ 0 • μ r • x • Tt . ( 11
)
The reluctances 23 and 32 are, respectively, given as

23 = 32 = Tch μ 0 • μ r • (Tw-x) 2 • Tt . ( 12
)

E. Reluctances of the Leakage Flux

The RN model developed in this paper has two advantages. Firstly, both flux leakages between the stator and the translator on one part and on the stator on the other are considered and represented by the circumference arcs. This leakage flux strongly influences the flux with a considerable effect on the magnetic field in this actuator. In fact, RN simulation results obtained by neglecting the leakage flux lead to important errors compared with the FEM. Secondly, the flux leakage reluctances vary with the translator displacement. The leakage flux is modeled by a magnetic reluctance given by the following expression:

f = θ μ 0 • Tt 1 ln R ex R in (13)
where θ is the angle of the circumferences arc. R ex and R in are, respectively, the exterior radius and interior radius of the circumferences arc, as shown in Fig. 4.

Each part of the flux leakage is modeled by an adapted magnetic reluctance. The parameters of the magnetic reluctances are parameterized with the displacement of the mobile part. Indeed, the flux path of each part increases or decreases with the displacement of the mobile part. In this RN model, the leakage flux reluctances change according the translator displacement. The RN is modeled in software Reluctool, and the finite element model of 1/4 LHM is shown in Fig. 5.

When the translator is moved from the unaligned position to the aligned position, the flux modeled using the reluctance f 2 decreased considerably. The angle formed is equal to 1.57 rad. The interior radius is constant and equal to the half of the air gap. However, the exterior radius decreases according the expression Twx. The f 2 is given by the following expression:

f 2 = 1, 57 μ 0 • Tt 1 ln Tw-x e/2 . ( 14
)
The exterior radius of the flux modeled using f 4 decreases according the expression (Twx)/2. Indeed, in the unaligned position, the exterior radius is equal to the half of the tooth width. Applying, the same methodology, all of flux leakage reluctances is determined. The f 3 and f 5 are, respectively, given as

f 3 = 1, 57 μ 0 • Tt 1 ln Tw e/2+x
(15)

f 5 = 1, 57 μ 0 • Tt 1 ln Tw/2 e/2+x/2 . ( 16
)
The f 4 and f 6 are written as

f 4 = f 6 = 1, 57 μ 0 • Tt 1 ln (Tw-x)/2 e/2 . ( 17
)
The reluctance f 1 represents the flux leakage canalized between the two poles in the same stator. The flux leakage in this case is constant. The angle formed is 3, 14 rad. Consequently, f 1 is given as

f 1 = 3, 14 μ 0 • Tt 1 ln 1,5•Tw e/2 . ( 18
)

III. STATIC FORCES COMPARISON

The algorithm that solves the RN is shown in Fig. 6. It allows reducing the computation time. A set of mobile part positions has to be defined before generating the model.

The propulsive forces calculated with RN model implemented by the Reluctool software and the FEM for three different current excitations are shown in Fig. 7. As shown earlier, force curves are almost identical. Hence, the RN model of the LHM is in accordance with the FEM.

IV. LHM OPTIMIZATION

The objective function of this optimization consists in minimizing the LHM mass while keeping the propulsive force. Moreover, the LHM mass is a very important parameter in the level of actuator manufacturing and transportation. However, many constraints must be respected. In particular, the propulsive force must be kept greater than 5000 N. The optimization constraints are summarized in Table II.

The RN model proposed in this paper reduces the computation time in both modeling and optimization of the LHM. The optimization process ends in 3463 ms after 28 iterations. The initial and optimal geometrical and electrical parameters are given in Table III.

The initial and optimized geometrical structures are shown in Fig. 8.

The mass of different parts of the initial and optimized structures is summarized in Table IV. The maximum propulsive force is also given in the same table.

The propulsive force of the optimized structure is calculated using two methods: 1) the RN model and 2) the FEM. The obtained results are shown in Fig. 9, which confirms the robustness of the developed RN model.

V. CONCLUSION

A precise RN model is developed to model, design, and optimize a linear actuator. This model is developed using the Reluctool software and is validated using the FEM. The efficiency of this model is its ability to generate automatically the electromagnetic characteristics for any mobile part positions. This model has a good numerical stability and leads to precise force computation. The model is applied for sizing by optimization. The sequential quadrating programming method is used, and accurate optimization results are obtained. The simulation results presented and discussed confirm the robustness of the developed RN model and the optimization approach.

Fig. 1 .

 1 Fig. 1. Double-sided hybrid motor.

Fig. 2 .

 2 Fig. 2. Geometrical dimensions.

Fig. 3 .

 3 Fig. 3. Construction of the RN.

Fig. 4 .

 4 Fig. 4. Geometrical view of leakage flux tube.

Fig. 5 .

 5 Fig. 5. RN model and finite element model of 1/4 LHM.

Fig. 6 .

 6 Fig. 6. RN model algorithm.

Fig.

  Fig. Initial and optimized LHM geometry.

TABLE I LHM DESIGN SPECIFICATIONS Efficient Reluctance Network Formulation for Modeling Design and Optimization of Linear Hybrid Motor

 I . M. Barhoumi 1,2,3 , F. Wurtz 3 , Member, IEEE, C. Chillet 3 , B. Ben Salah 1 , and O. Chadebec 3

EI. INTRODUCTION

TABLE II OPTIMIZATION CONSTRAINTS TABLE III OPTIMIZATION

 IICONSTRAINTSIII RESULTS OF GEOMETRICAL PARAMETERS

TABLE IV OPTIMIZATION

 IV RESULTS OF MASS AND PROPULSIVE FORCE9. Static force of the optimal structure calculated with FEM and RN.

ACKNOWLEDGMENT

This work is supported by University of Tunis El Manar (Tunisia), Univ. Grenoble Alpes, G2Elab (France), and Almajmaah University (Saudi Arabia).