J.-P Guillement
email: guillement@math.univ-nantes.fr

R G Novikov
email: novikov@cmap.polytechnique.fr

Inversion of weighted Radon transforms via finite Fourier series weight approximations

We consider weighted Radon transforms on the plane. We show that the Chang approximate inversion formula for these transforms admits a principal refinement as inversion via finite Fourier series weight approximations. We illustrate this inversion approach by numerical examples for the case of the attenuated Radon transforms in the framework of the single-photon emission computed tomography (SPECT).

Introduction

A basic problem of many tomographies consists in finding an unknown function f on R 2 from its weighted ray transform P W f on R × S 1 for some known weight W , where

P W f (s, θ) = R W (sθ ⊥ + tθ, θ)f (sθ ⊥ + tθ)dt, s ∈ R, θ = (θ 1 , θ 2 ) ∈ S 1 , θ ⊥ = (-θ 2 , θ 1 ), (1.1)
where f = f (x), W = W (x, θ), x ∈ R 2 . Up to change of variables, the operator P W is known also in the literature as the weighted Radon transform on the plane. In this work we always assume that

W ∈ C(R 2 × S 1 ) ∩ L ∞ (R 2 × S 1 ), w 0 (x) def = 1 2π S 1 W (x, θ)dθ = 0, x ∈ R 2 , (1.2)
where W is complex-valued, in general, and dθ is the standard element of arc length on S 1 . Additional assumptions on W will be formulated in the framework of context. In particular, in important particular cases W is real-valued and strictly positive:

W = W , W ≥ c > 0.
(1.3)

In addition, for this work one can always assume that

f ∈ L ∞ (R 2 ), supp f ⊂ D, (1.4)
where f is complex-valued, in general, and D is an open bounded domain (which is fixed a priori).

In definition (1.1) the product R × S 1 is interpreted as the set of all oriented straight lines in R 2 . If γ = (s, θ) ∈ R × S 1 , then γ = {x ∈ R 2 : x = sθ ⊥ + tθ, t ∈ R} (modulo orientation) and θ gives the orientation of γ.

If W ≡ 1, then P W is reduced to the classical ray (or Radon) transform on the plane. If W (x, θ) = exp -+∞ 0 a(x + tθ)dt ,

(1.5)

where a is a complex-valued sufficiently regular function on R 2 with sufficient decay at infinity, then P W is known in the literature as the attenuated ray (or Radon) transform on the plane. The classical ray transform arises, in particular, in X-ray transmission tomography, see e.g. [Na]. The attenuated ray transform (at least, with real-valued a ≥ 0) arises, in particular, in single-photon emission computed tomography (SPECT), see e.g. [Na], [Ku]. Transforms P W with some other weights also arise in applications, see e.g. [Q], [Ku], [MP], [Ba].

Exact and simultaneously explicit inversion formulas for the classical and attenuated ray (or Radon) transforms on the plane were given for the first time in [Ra] and [START_REF] Novikov | An inversion formula for the attenuated X-ray transformation[END_REF], respectively. For some other weights W , exact and simultaneously explicit inversion formulas were given in [TM], [BS], [G], [START_REF] Novikov | Weighted Radon transforms for which Chang's approximate inversion formula is exact[END_REF]. Note that we say that an inversion method for P W is an explicit inversion formula if its complexity is comparable with complexity of the aforementioned Radon inversion formula of [Ra].

Apparently, for general P W , under assumptions (1.2), (1.3), explicit and simultaneously exact (modulo Ker P W ) inversion formulas are impossible. In addition, [Bo] gives an example of infinitely smooth W satisfying (1.2), (1.3) and such that Ker P W = 0 in the space of infinitely smooth compactly supported functions on R 2 . But, due to [BQ], Ker P W = 0 in the space of continuous compactly supported functions on R 2 for realanalytic W satisfying (1.3).

On the other hand, the following Chang approximate (but explicit) inversion formula for P W , where W is given by (1.5) with a ≥ 0, has been used for a long time (see e.g. [Ch], [Ku], [START_REF] Novikov | Weighted Radon transforms for which Chang's approximate inversion formula is exact[END_REF]):

f appr = (w 0 ) -1 P -1 P W f, (1.6)
where w 0 is defined in (1.2), P -1 denotes the classical Radon inversion realized via the formula

P -1 q(x) = 1 4π S 1 θ ⊥ ∇ x 1 π p.v. R q(t, θ) xθ ⊥ -t dt dθ, x ∈ R 2 , θ ⊥ ∇ x = -θ 2 ∂/∂x 1 + θ 1 ∂/∂x 2 , xθ ⊥ = -θ 2 x 1 + θ 1 x 2 , θ = (θ 1 , θ 2 ), x = (x 1 , x 2 ), (1.7)
where q is a test function on R × S 1 . It is known that formula (1.6) is efficient for the first approximation in SPECT reconstructions and, in particular, is sufficiently stable for reconstructions from discrete data with strong Poisson noise. The exact inversion formula of [START_REF] Novikov | An inversion formula for the attenuated X-ray transformation[END_REF] for the attenuated ray transform is considerably less stable in this respect. For more information in this connection see [START_REF] Guillement | Optimized analytic reconstruction for SPECT[END_REF]. Formula (1.6), under assumptions (1.2), can be considered as approximate inversion of P W via the approximation (1.8) i.e. via the zero term Fourier approximation for W in angle variable. In addition, due to [Ra], formula (1.5) is exact if

W (x, θ) ≈ w 0 (x), x ∈ R 2 , θ ∈ S 1 ,
W (x, θ) ≡ w 0 (x), x ∈ R 2 , θ ∈ S 1 .
(1.9)

Besides, due to [START_REF] Novikov | Weighted Radon transforms for which Chang's approximate inversion formula is exact[END_REF], formula (1.5) is exact if and only if

W (x, θ) -w 0 (x) ≡ w 0 (x) -W (x, -θ), x ∈ R 2 , θ ∈ S 1 .
(1.10) However, already for W of (1.5) property (1.10) is not fulfilled, in general.

In the present work we consider approximate inversion of P W via the approximations

W (x, θ(ϕ)) ≈ N n=-N e inϕ w n (x), w n (x) = 1 2π π -π e -inϕ W (x, θ(ϕ))dϕ, x ∈ R 2 , θ(ϕ) = (cos ϕ, sin ϕ), ϕ ∈ [-π, π], N ∈ N ∪ 0.
(1.11)

One can see that for N = 0 approximation (1.11) is reduced to (1.8). Our approximate inversion algorithms for P W on functions of (1.4), under assumptions (1.2), are presented in Subsection 2.4 of Section 2. In these considerations we proceed from [Ku] and [START_REF] Novikov | Weighted Radon transforms and first order differential systems on the plane[END_REF].

In Section 3, our approximate inversions (2.20), (2.22) of Subsection 2.4 are illustrated by numerical examples for W given by (1.5) with a ≥ 0 in the framework of SPECT.

Approximate inversion of P

W 2.1. Symmetrization of W . Let A W f = P -1 P W f, (2.1)
where P -1 is defined by (1.7). Let

W sym (x, θ) = 1 2 (W (x, θ) + W (x, -θ)), x ∈ R 2 , θ ∈ S 1 . (2.2)
We consider also the Fourier series

W (x, θ(ϕ)) = +∞ n=-∞ e inϕ w n (x), x ∈ R 2 , ϕ ∈ [-π, π], (2.3)
where w n are defined in (1.11), θ(ϕ) = (cos ϕ, sin ϕ).

The following formulas hold (see [Ku], [START_REF] Novikov | Weighted Radon transforms for which Chang's approximate inversion formula is exact[END_REF]):

1 2 P W f (s, θ) + P W f (-s, -θ) = P W sym f (s, θ), (s, θ) ∈ R × S 1 , (2.4) A W f = P -1 P W sym f, (2.5) W sym (x, θ(ϕ)) = +∞ l=-∞ e i2lϕ w 2l (x), x ∈ R 2 , ϕ ∈ [-π, π].
(2.6)

Actually, using (2.4)-(2.6) we reduce inversion of P W to inversion of P W sym . In particular, using (2.1), (2.5), (2.6), one can see that such a reduction arises already in the framework of (1.6).

Operators

Q W,D,m and numbers σ W,D,m . Consider z = x 1 + ix 2 , z = x 1 -ix 2 , x = (x 1 , x 2 ) ∈ R 2 .
(2.7)

Using (2.7) we identify R 2 with C. Let Π, Π denote the linear integral operators on R 2 identified with C such that Πu(z) = - 1 π C u(ζ) (ζ -z) 2 dRe ζ dIm ζ, Πu(z) = - 1 π C u(ζ) ( ζ -z) 2 dRe ζ dIm ζ, (2.8)
where u is a test function, z ∈ C; see e.g. [V] for detailed properties of these operators.

Let D be the domain of (1.4). Let χ D denote the characteristic function of D, i.e.

χ D ≡ 1 on D, χ D ≡ 0 on R 2 \D. (2.9) Let Q W,D,m = m l=1 (-Π) l w 2l w 0 + (-Π) l w -2l w 0 χ D for m ∈ N, Q W,D,m = 0 for m = 0, (2.10)
where w 0 , w 2l are the Fourier coefficients of (2.3), (2.6) and w ±2l /w 0 , χ D are considered as the multiplication operators on R 2 . Let

σ W,D,m = m l=1 sup x∈D w 2l (x) w 0 (x) + sup x∈D w -2l (x) w 0 (x) for m ∈ N, σ W,D,m = 0 for m = 0.
(2.11)

According to [START_REF] Novikov | Weighted Radon transforms and first order differential systems on the plane[END_REF] we have that

Q W,D,m L 2 (R 2 )→L 2 (R 2 ) ≤ σ W,D,m .
(2.12) 2.3. Exact inversion for finite Fourier series weights. Let conditions (1.2), (1.4) be fulfilled and

W (x, θ(ϕ) = N n=-N e inϕ w n (x), x ∈ R 2 , θ(ϕ) = (cos ϕ, sin ϕ), ϕ ∈ [-π, π], N ∈ N ∪ 0. (2.13) Suppose also that σ W,D,m < 1 for m = [N/2], (2.14)
where [N/2] denotes the integer part of N/2. Then according to [START_REF] Novikov | Weighted Radon transforms and first order differential systems on the plane[END_REF] we have the following exact inversion for P W :

f = (w 0 ) -1 (I + Q W,D,m ) -1 P -1 P W f, (2.15)
where I is the identity operator on R 2 , P -1 is defined by (1.7). In addition, in view of (2.12), (2.14) we have that

(I + Q W,D,m ) -1 = I + +∞ j=1 (-Q W,D,m ) j .
(2.16)

In addition, under our assumptions, we have that (see [START_REF] Novikov | Weighted Radon transforms and first order differential systems on the plane[END_REF]):

A W f = P -1 P W f ∈ L 2 (R 2 ).
(2.17) 2.4. Approximate inversion for general weights. Let conditions (1.2), (1.4) be fulfilled and

∞ l=1 w 2l w 0 L 2 (D) + w -2l w 0 L 2 (D) < +∞, (2.18)
where w n are defined by (2.11). Suppose that

σ W,D,m < 1 for fixed m ∈ N ∪ 0. (2.19)
Then we consider that

f ≈ f m , (2.20a) f m = (w 0 ) -1 (I + Q W,D,m ) -1 P -1 P W f, (2.20b)
where I is the identity operator, Q W,D,m is defined by (2.10), P -1 is defined by (1.7). Note that the right-hand side of (2.20b) is well-defined: in particular, in view of (2.12), (2.19) we have formula (2.16) in norm • L 2 (R 2 )→L 2 (R 2 ) and using (2.18) one can show that

P -1 P W f ∈ L 2 (R 2 ).
(2.21) Formula (2.20) is a natural extension of the Chang formula (1.6). In particular, (2.20) for m = 0 coincides with (1.6).

In addition, if (2.19) is fulfilled for some m ≥ 1, then f m is a refinement of the Chang approximation f 0 and, more generally, f j is a refinement of

f i for 0 ≤ i < j ≤ m. But of course f j = f i if w 2l ≡ 0, w -2l ≡ 0 for i < l ≤ j.
Actually, in many important examples condition (2.19) is efficiently fulfilled for small m (e.g. m = 2) and is not fulfilled for great m. Therefore, in the present work we propose the following approximate reconstruction of f from P W f : find maximal m such that (2.19) is still efficiently fulfilled and approximately reconstruct f via (2.20) using (2.12), (2.16).

(2.22)

In Section 3 we illustrate this approximate reconstruction by numerical examples in the framework of SPECT. In particular, these numerical examples show that stability properties of this approximate reconstruction f m are similar to stability properties of the Chang approximate reconstruction f 0 but f m is more precise than f 0 if m > 0.

Using considerations of Subsection 2.3 one can see also that

f = f m -(w 0 ) -1 (I + Q W,D,m ) -1 P -1 P δW m f, (2.23)
where

δW m (x, θ(ϕ)) def = W (x, θ(ϕ)) - 2m+1 n=-2m-1 e inϕ w n (x), (2.23) x ∈ R 2 , θ(ϕ) = (cos ϕ, sin ϕ), ϕ ∈ [-π, π], m ∈ N ∪ 0, (2.24)
where w n are defined by (1.11). One can use (2.23) for estimating the error f -f m . One can also consider (2.23) as an integral equation for finding f from f m . Note also that equation (2.23) for m = 0 is, actually, well-known, see e.g. [Ku].

2.5. Relations with [Ku]. Let conditions (1.2), (1.4) be fulfilled and

lim m→+∞ σ W,D,m = σ W,D < +∞.
(2.25)

Then we can consider

Q W,D = lim m→+∞ Q W,D,m in • L 2 (R 2 )→L 2 (R 2 ) .
(2.26)

In addition, if σ W,D < 1, (2.27) then f = (w 0 ) -1 (I + Q W,D ) -1 P -1 P W f, (2.28) (I + Q W,D ) -1 = I + +∞ j=1 (-Q W,D ) j .
(2.29)

Actually, (2.28) is a linear integral equation for exact reconstruction of f from P W f under assumptions (1.2), (1.4), (2.27). In addition, (2.29) can be interpreted as the method of successive approximations for solving this equation. It is possible to show that the reconstruction of f from P W f of [Ku] (or, more precisely, the linear integral equation on ε on page 814 of [Ku]) can be transformed into (2.28).

In [Ku] (under the assumption that 0 < W < 1) it was shown that this equation on ε of [Ku] is uniquely solvable if

lim m→+∞ ρ W,D,m = ρ W,D < 1, ρ W,D,m def = m l=1 (sup x∈D |w 2l (x)| + sup x∈D |w -2l (x)|) min x∈D |w 0 (x)| .
(2.30)

One can see that in many cases σ W,D,m is much smaller than ρ W,D,m and condition (2.27) is much less restrictive than (2.30). Actually, in order to relate considerations of [Ku] on one hand with considerations of [START_REF] Novikov | Weighted Radon transforms and first order differential systems on the plane[END_REF] and the present work on the other hand, we use the following formulas

F(Π l u)(ξ) = ξ 1 -iξ 2 ξ 1 + iξ 2 l Fu(ξ), F( Πl u)(ξ) = ξ 1 + iξ 2 ξ 1 -iξ 2 l Fu(ξ), ξ = (ξ 1 , ξ 2 ) ∈ R 2 , l ∈ N, (2.31) Q W,D,m u(x) = P -1 (P W,D,m u)(x), P W,D,m u(s, θ) def = R -1 l=-m + m l=1 w 2l (sθ ⊥ + tθ) w 0 (sθ ⊥ + tθ) (θ 1 + iθ 2 ) 2l × χ D (sθ ⊥ + tθ)u(sθ ⊥ + tθ)dt, x ∈ R 2 , s ∈ R, θ = (θ 1 , θ 2 ) ∈ S 1 , m ∈ N, (2.32)
where F denotes the 2D Fourier transform operator, Π, Π are defined by (2.8), Q m is defined by (2.10), P -1 is defined by (1.7), u is a test function.

Finally, note that in our numerical examples of Section 3 the approximate reconstruction (2.20) is realized numerically using formula (2.32).

Numerical examples

Framework of SPECT.

All numerical examples of this work are given in the framework of SPECT. In particular, we assume that W is given by (1.5), where a ≥ 0. We recall that in SPECT, after restricting the problem to a fixed 2D plane, the functions P W f , a, f are interpreted as follows:

• f is distribution of radioactive isotopes emitting photons;

• a is photon attenuation coefficient;

• in addition, it is assumed that supp a ⊆ D, supp f ⊆ D, where D is some known compact domain;

• P W f describes the expected radiation outside D.

In addition, in some approximation, measured SPECT data are modeled as P W f with Poisson noise on some discrete subset Γ of R × S 1 . Usually, it is also assumed that:

• D ⊂ B R = {x ∈ R 2 : |x| ≤ R}, where R is the radius of image support, • Γ is a uniform n × n sampling of T R = {γ ∈ R × S 1 : γ ∩ B R = 0} = {(s, θ) ∈ R × S 1 : |s| ≤ R}. (3.1)
For more information on SPECT see e.g. [Na], [Br], [START_REF] Guillement | On Wiener type filters in SPECT[END_REF] and references therein.

3.2. Image parameters. All 2D images of this work are considered on n × n grids (X and Γ), where n = 128. We assume that X is a uniform n × n sampling of

D R = {x = (x 1 , x 2 ) ∈ R 2 : |x 1 | ≤ R, |x 2 | ≤ R} (3.2)
and Γ is a uniform n × n sampling of T R , where R is the radius of the image support.

In addition, all 2D images of this work are drawn using a linear greyscale in such a way that the dark grey color represents zero (or negative values, if any) and white corresponds to the maximum value of the imaged function.

3.3. Numerical phantoms. We consider two numerical phantoms (phantom 1 and phantom 2). Attenuation maps a, emitter activities f , noiseless emission data g = P W f and noisy emission data p for these phantoms are shown in figures 1 and 2. Phantom 1 is a version of the elliptical chest phantom (used for numerical simulations of cardiac SPECT imaging; see e.g. [Br], [START_REF] Guillement | On Wiener type filters in SPECT[END_REF], [START_REF] Guillement | Optimized analytic reconstruction for SPECT[END_REF]). Actually, this version is the same as in [START_REF] Guillement | On Wiener type filters in SPECT[END_REF], [START_REF] Guillement | Optimized analytic reconstruction for SPECT[END_REF] and, in addition to figure 1, its description includes the following information:

• the major axis of the ellipse representing the body is 30 cm, • attenuation coefficient a is 0.04 cm -1 in the lung region (modeled as two interior ellipses), 0.15 cm -1 elsewhere within the body ellipse, and zero outside the body,

• emitter activity f is in the ratio 8:0:1:0 in myocardium (represented as a ring), lungs, elsewhere within the body, and outside the body, • noisy emission data p contain 30 percents of the Poisson noise in the sense of L 2 -norm and the total number of photons is 125 450.

Phantom 2 is a simulated numerical version of the Utah phantom (designed at the 2 nd International Meeting on fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, Snowbird, Utah, 1993). A real non simulated version of this phantom was considered, in particular, in [GJKNT]. However, in the present work we consider its simulated numerical version in order to see clearly different reconstruction errors. In addition to figure 2, the description of phantom 2 includes the following information:

• geometrically the phantom consists of a large disk containing two small disks, where the radius of the large disks 10 cm,

• attenuation coefficient a is 0.16 cm -1 in the large disk outside the small disks, 0.63 cm -1 in the left small disk, 0.31 cm -1 in the right small disk, and zero outside the large disk,

• emitter activity f is a positive constant in the large disk outside the small disks and zero elsewhere • noisy emission data p contain 30 percents of the Poisson noise in the sense of L 2 -norm and the total number of photons is 89 350.

Note also that all computations of the present work are fulfilled using the same numerical realizations of basic formulas as in [GJKNT], [START_REF] Guillement | On Wiener type filters in SPECT[END_REF], [START_REF] Guillement | Optimized analytic reconstruction for SPECT[END_REF].

m = 1| m = 2| m = 3| m = 20| σ W,D,m | 0.489 | 0.694 | 0.803 | 1.296 | ρ W,D,m | 3.112 | 4.567 | 5.436 | 9.719 |
For phantoms 1 and 2, tables 1 and 2 show that condition (2.19) is efficiently fulfilled, at least, for m = 1 and m = 2, whereas ρ W,D,m > 1 already for m = 1.

Reconstruction results.

For phantoms 1 and 2 we consider the approximate reconstructions f m realized numerically via (2.20b) (using the method of successive approximations i.e. using (2.16)) from the noiseless data g and filtered noisy data p. In addition:

• we consider f m , mainly, under the condition that σ W,D,m ≤ 0.7,

• we consider p = Wp for W = A sym 8,8 , where W = A sym l 1 ,l 2 is the approximately optimal space-variant Wiener-type filter of Section 5.3 of [START_REF] Guillement | On Wiener type filters in SPECT[END_REF].

In addition to f m we consider also their non-negative parts f + m , where

f + m (x) = f m (x) if f m (x) ≥ 0 and f + m (x) = 0 if f m (x) < 0
For phantoms 1 and 2, tables 1 and 2 show that σ W,D,m ≤ 0.7 is fulfilled for m = 0, 1, 2 only.

To study reconstruction errors we consider, in particular,

η(u 2 , u 1 , X) = u 2 -u 1 L 2 (X) u 1 L 2 (X) , (3.3)
where u 1 , u 2 are test functions on X.

The approximate reconstructions f 0 , f 2 with their central horizontal profiles from the noiseless data g for phantoms 1 and 2 are shown in figures 3 and 4.

Figures 3,4 and tables 3, 4 show that for phantoms 1 and 2 for the noiseless case the approximations f 1 , f 2 are considerably more precise than the classical Chang approximation f 0 .

The approximate reconstructions f 0 , f 2 with their central horizontal profiles from the filtered noisy data p = A sym 8,8 p for phantoms 1 and 2 are shown in figures 5 and 6, where A sym l 1 ,l 2 is the aforementioned filter of [START_REF] Guillement | On Wiener type filters in SPECT[END_REF]. Tables 5 and6 show the relative errors η(f m , f, X) and η(f + m , f, X) in L 2 -norm for f m reconstructed from filtered noisy data p = A sym 8,8 p for phantoms 1 and 2. Table 5. Relative errors η(f m , f, X), η(f + m , f, X) for f m reconstructed from p = A sym 8,8 p for phantom 1.

Figures 5, 6 and tables 5, 6 show that for phantoms 1 and 2 for the noisy case the approximations f 1 , f 2 are also more correct than the classical Chang approximation f 0 .

Conclusions

In this work we showed that the classical Chang approximate inversion formula (1.6) admits a very natural extention into inversion via finite Fourier series weight approximations or, more precisely, into inversion via (2.20) considered under assumption (2.19). Related theoretical considerations are presented in Sections 1 and 2 and numerical examples in the framework of SPECT are given in Section 3. Our examples of Section 3 include comparisons with the approximate Chang reconstruction f 0 and show numerical efficiency (with respect to precision and stability) of our approximate reconstructions f m for m > 0, under the condition that inequality (2.19) is efficiently fulfilled.

Note also that considerations of Subsections 2.5 and 3.4 explain convergence of the iterative inversion of [Ku] for many cases when the inequality of (2.30) is not fulfilled. The point is that less restrictive inequality (2.27) is sufficient for this convergence.

Figure 1 .Figure 2 .

 12 Attenuation map a (a), emitter activity f (b), noiseless emission data g = P W f (c), noisy emission data p (d) for phantom 1. (See Subsections 3.1, 3Attenuation map a (a), emitter activity f (b), noiseless emission data g = P W f (c), noisy emission data p (d) for phantom 2. (See Subsections 3.1, 3.3)

3. 4 .

 4 Results for the bounds σ W,D,m and ρ W,D,m . The bound numbers σ W,D,m , ρ W,D,m of (2.11), (2.30) for phantoms 1 and 2 are shown in tables 1 and 2, where D = B R .

Figure 5 .Figure 6 .

 56 Approximate reconstructions f 0 (a) and f 2 (c) with their central horizontal profiles (b) and (d) from the filtered noisy data p = A sym 8,8 p for phantom 1. Approximate reconstructions f 0 (a) and f 2 (c) with their central horizontal profiles (b) and (d) from the filtered noisy data p = A sym 8,8 p for phantom 2.

Table 1 .

 1 Numbers σ W,D,m , ρ W,D,m of (2.11), (2.30) for phantom 1, where D = B R .

		m = 1|	m = 2|	m = 3|	m = 20|
	σ W,D,m |	0.390 |	0.584 |	0.739 |	1.320 |
	ρ W,D,m |	1.399 |	2.025 |	2.494 |	4.532 |

Table 2 .

 2 Numbers σ W,D,m , ρ W,D,m of (2.11), (2.30) for phantom 2, where D = B R .

Acknowledgements

The second author was partially supported by TFP No 14.A18.21.0866 of Ministry of Education and Science of Russian Federation (at Faculty of Control and Applied Mathematics of Moscow Institute of Physics and Technology).

 3 and4show the relative errors η(f m , f, X) and η(f + m , f, X) in L 2 -norm for the approximate reconstructions f m with respect to precise f for the noiseless case for phantoms 1 and 2.

Table 3. Relative reconstruction errors η(f m , f, X), η(f + m , f, X) for the noiseless case for phantom 1.


Simulating the effects of spatial configurations of agricultural ditch drainage networks on surface runoff from agricultural catchments

Keywords: Agricultural ditch drainage network, cultivated catchments, stochastic simulations, hydrological modeling, runoff

The study of runoff is a crucial issue because it is closely related to flooding, water quality and erosion. In cultivated catchments, agricultural ditch drainage networks are known to influence runoff. As anthropogenic elements, agricultural ditch drainage networks can therefore be altered to better manage surface runoff in cultivated catchments. However, the relationship between the spatial configuration, i.e., the density and the topology, of agricultural ditch drainage networks and surface runoff in cultivated catchments is not understood. We studied this relationship by using a random network simulator that was coupled to a distributed hydrological model. The simulations explored a large variety of spatial configurations corresponding to a thousand stochastic agricultural ditch drainage networks on a 6.4 km² Mediterranean cultivated catchment. Next, several distributed hydrological functions were used to compute water flow-paths and runoff for each simulation.

The results showed that (i) denser networks increased the drained volume and the peak discharge and decreased hillslopes runoff, (ii) greater network density did not affect the surface runoff any further above a given network density, (iii) the correlation between network density and runoff was weaker for small subcatchments (< 2 km²) where the variability in the drained area that resulted from changes in agricultural ditch drainage

Introduction

Hedges, ditches, terrace fronts, embankments, grass strips and roads are common linear features of cultivated landscapes (Figure 1).

Figure 1

The functional importance of these features has been emphasized in various contexts and biophysical processes. In ecology, they provide valuable habitats, enhance connectivity within landscapes and can serve as habitats for rare species [START_REF] Forman | Hedgerows and hedgerow networks in landscape ecology[END_REF][START_REF] Herzon | Agricultural drainage ditches, their biological importance and functioning[END_REF][START_REF] Watson | The distribution of three uncommon freshwater gastropods in the drainage ditches of British grazing marshes[END_REF][START_REF] Pita | Conserving the Cabrera vole, Microtus cabrerae, in intensively used Mediterranean landscapes[END_REF]. Linear features affect soil redistribution as well, especially terraces that attenuate slopes, facilitate agriculture and limit long-term soil loss [START_REF] Bevan | Terraced fields and Mediterranean landscape structure: An analytical case study from Antikythera, Greece[END_REF][START_REF] Paroissien | A regional-scale study of multidecennial erosion of vineyard fields using vine-stock unearthing-burying measurements[END_REF]. By favoring water infiltration, they play a role in groundwater hydrology (Dages et al., 2009). Linear features also alter overland flow paths [START_REF] Duke | Incorporating ancillary data to refine anthropogenically modified overland flow paths[END_REF][START_REF] Gascuel-Odoux | Incorporating landscape features to obtain an object-oriented landscape drainage network representing the connectivity of surface flow pathways over rural catchments[END_REF]. They can force water to not follow the steepest slope, but they can also concentrate runoff along the steepest slope.

When focusing on the impacts of agricultural ditch drainage networks, four roles are commonly considered: the interception of overland flow on hillslopes, the drainage of groundwater and exfiltration to ditch networks by lowering the water table, infiltration from the ditch towards the groundwater and conveyance of water towards downstream areas [START_REF] Adamiade | Influence d`un fossé sur les écoulements rapides au sein d`un versant[END_REF][START_REF] Carluer | Assessment and modelling of the influence of man-made networks on the hydrology of a small watershed: implications for fast flow components, water quality and landscape management[END_REF][START_REF] Dunn | Modelling the hydrological impacts of open ditch drainage[END_REF]. A role will be predominant or will not, depending on various factors, including climate, soil, direction to the steepest slope and the location of the network within a catchment.

In Mediterranean catchments, where short but intense storms predominate, Hortonian runoff is much more important than subsurface flow [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF]. Hence, considering surface runoff is crucial because it is closely related to flooding, water quality and erosion [START_REF] Fiener | Spatio-temporal patterns in land use and management affecting surface runoff response of agricultural catchments--A review[END_REF].

To control runoff in agricultural catchments, we can work with many anthropogenic elements, such as land use and the tillage practices (Colin et al., 2011 a;[START_REF] Souchère | Modelling the impact of agri-environmental scenarios on runoff in a cultivated catchment (Normandy, France)[END_REF][START_REF] Takken | The prediction of runoff flow directions on tilled fields[END_REF] or the use of grassed waterways [START_REF] Fiener | Effectiveness of grassed waterways in reducing runoff and sediment delivery from agricultural watersheds[END_REF]. Agricultural ditch drainage networks are also one of the few structural landscape elements that control runoff and that can be changed without consuming or drastically modifying the agricultural area. Therefore, the relationship between agricultural ditch drainage networks and runoff is a crucial landscape structure-function question that is worth investigating. In addition to this relationship, the optimization of agricultural ditch drainage networks is of interest to agricultural landscape managers, for example, in minimizing floods or the fate of pesticides (DGFAR, 2008).

Consequently, such an investigation of the relationship between agricultural ditch drainage network spatial configurations and runoff will permit an assessment of the hydrological benefits that we could expect through modifying these networks.

Because of the emphasized role of ditches in cultivated landscape hydrology, an increasing number of spatially distributed models have explicitly considered ditch networks and their spatial configurations in landscapes, which have always been considered to be directed tree structures [START_REF] Dha | Hydrological modelling of a drained grazing marsh under agricultural land use and the simulation of restoration management scenarios[END_REF][START_REF] Branger | Utilisation d'une plate-forme de modélisation environnementale pour représenter le rôle d'aménagements hydro-agricoles sur les flux d'eau et de pesticides. Application au bassin versant de la Fontaine du Theil[END_REF][START_REF] Carluer | Assessment and modelling of the influence of man-made networks on the hydrology of a small watershed: implications for fast flow components, water quality and landscape management[END_REF][START_REF] Dunn | Modelling the hydrological impacts of open ditch drainage[END_REF][START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF]. Consequently, these models allow for the assessment of the hydrological impact of agricultural ditch drainage networks. These models can also be used to test certain scenarios concerning possible evolutions of these networks. For example, [START_REF] Krause | Modelling the impacts of land-use and drainage density on the water balance of a lowland-floodplain landscape in northeast Germany[END_REF] modeled the effect of removing part of an agricultural ditch drainage network on the hydrology of a lowland floodplain in northeast Germany and showed that groundwater recharge was altered. In the Mediterranean area, [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF] found that the substitution of a ditch network by a hypothetical natural drainage network that followed the steepest slopes could affect runoff by decreasing peak discharge and increasing lag time.

However, these latter studies only considered a few scenarios, such as the suppression of the agricultural parts of the drainage network.

Because these studies only considered a limited variability of the ditch network, they could not determine, even for a given catchment, the relationship between the agricultural ditch drainage network spatial configuration, i.e., network density and topology, a landscape structural property, and surface runoff, a landscape functioning property. To explore the relationship between the structure and the function of a landscape, the usual way to proceed in environmental sciences is to couple stochastic landscape structure simulations to a landscape functioning model [START_REF] Colin | Small catchment agricultural management using decision variables defined at catchment scale and fuzzy rule based system: A Mediterranean vineyard case study[END_REF][START_REF] Gumiere | Contribution à la modélisation déterministe spatialisée de l'érosion hydrique des sols à l'échelle des petits bassins versants cultivés[END_REF][START_REF] Ber | Neutral modelling of agricultural landscapes by tessellation methods -Application for gene flow simulation[END_REF][START_REF] Van Nieuwenhuyse | Pattern-process relationships in surface hydrology: hydrological connectivity expressed in landscape metrics[END_REF][START_REF] Viaud | Modeling the impact of the spatial structure of a hedge network on the hydrology of a small catchment in a temperate climate[END_REF]. To date, there has not been an equivalent study of the impacts of agricultural ditch drainage networks on runoff.

In this study, the objective was to analyze the extent to which the spatial configuration and especially the density of an agricultural ditch drainage network could control the surface runoff of a given cultivated catchment. The investigated case study was the 6.4 km² Bourdic catchment, which is located in the Languedoc vineyards in southern France, where previous hydrological modeling studies could be used as source of data for the hydrological parameterization. We performed this investigation by coupling a ditch network simulator [START_REF] Bailly | A spatial stochastic algorithm to reconstruct artificial drainage networks from incomplete network delineations[END_REF] to the surface runoff functions of the MHYDAS model. We focused on the hydrological response both at the catchment and subcatchment outlets and over hillslopes to investigate the influence of the agricultural ditch drainage network spatial configuration, i.e., density and topology, on surface runoff at various spatial scales.

Material and Methods

Study area

The study area was the 6.4 km² Mediterranean Bourdic catchment (Figure 2). It is located 50 km west of Montpellier in southern France. The altitude varies between 45 m at the outlet to 128 m westwards. The southern part of the catchment is compounded by rather well delineated small subcatchments with gentle hillslopes whereas the northern part is relatively flat. The land cover consists of mainly vineyards, with a small amount of cereal fields and shrubs. The actual agricultural ditch drainage network is 72 km long and covers all of the catchment except for the limestone uplands (cuesta) located in the center of the catchment. In addition to the drainage networks that are mapped in hydrographic databases (BD TOPO ®, BD CARTHAGE ®), the ditches highly extend past the drainage network (72 km by 10.8 km). The actual ditch drainage network was surveyed during the summer 2010 at an average rate of 1.5 to 3 km² per day per person, depending on the difficulty of the terrain. Fiftycentimeter resolution aerial photographs from IGN (Insitut Géographique National) were used to locate the ditches. Agricultural ditch drainage networks are rarely represented in hydrographic databases. An exhaustive survey of the ditch sizes in a sub-catchment of 1 km² revealed that 75 % of the ditches had an upper width between 50 and 120 cm and a depth between 30 and 80 cm.

The Bourdic catchment contains the Roujan subcatchment (approximately 1 km²), where numerous studies have been conducted for twenty years on soil and water resources. Because it is in the Mediterranean area, high-intensity and short-duration storms are frequent and Hortonian overland flow dominates subsurface flow [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF]. This latter fact motivated the choice of this catchment for studying the impact of agricultural ditch drainage networks on runoff.

Summary of the methods

The general methodology of the paper relied on the coupling of a landscape simulator with a distributed hydrological model (Figure 3).

Figure 3

Stochastic simulations of ditch drainage networks

The ditch network simulator has previsously been described in detail in [START_REF] Bailly | A spatial stochastic algorithm to reconstruct artificial drainage networks from incomplete network delineations[END_REF], and only the main principles are described below. It uses as its inputs the lattice of the field units' boundaries and a Digital Terrain Model (DTM). Each segment of this lattice determines a potential location for a ditch. The role of the algorithm is to select whether each potential location is a ditch or not according to a stochastic drainage process and DTM uncertainties.

Start and end nodes of each segment of the lattice have an elevation value that is extracted from the DTM. This lattice is partially directed by providing a unique direction to the lattice segments that have start node and end nodes with significant differences in elevation, i.e., that exceed a parameter of elevation uncertainty. For these segments, the slope is calculated as the difference between the start node and the end node altitudes, divided by the length of the segment. If the difference in the nodes' altitudes does not exceed the elevation uncertainty, the direction of the segment is not fixed and may be different for different simulations. The segments of the lattice that correspond to the main downstream parts of the network are considered invariant and are simulated as ditches for all of the ditch network simulations.

They define the invariant main downstream network.

The method of network generation consists of a stochastic drainage-like algorithm. It generates directed tree network structures corresponding to connected sub-graphs of the directed lattice of the agricultural field units' boundaries. The method is based on (i) directed random walks throughout the directed lattice of the field units that connect randomly selected segments to the invariant main downstream segments and (ii) a random branching/pruning process enabling the convergence to a targeted network length. The only parameters of this algorithm include the target network length, the tolerance parameter and the elevation uncertainty parameter. Once the simulation process leads to a simulated network with a network length that is equal to the target network length greater than or less than the tolerance parameter, the simulation process stops and the simulated network is saved. With this stochastic simulation process, each field unit boundary can be a ditch in a given simulation.

Numerous networks can thus be simulated to represent a wide variability of spatial configurations, i.e., density and topology.

In our case study, the invariant main downstream segments came from the French national databases on hydrography (BD TOPO ®, BD CARTHAGE ®), representing in total 10.8 km of segments (Figure 2). In fact, these segments mainly corresponded to actual channelized streams and not to ditches. However, for purposes of simplification, we called agricultural ditch drainage networks (or ditch networks) our simulated networks even if they contained these streams. A 5 m resolution photogrammetric DTM 1 was used for the network simulation process. The elevation uncertainty parameter was fixed to 1 m according to the accuracy of the DTM. In contrast to the study of [START_REF] Bailly | A spatial stochastic algorithm to reconstruct artificial drainage networks from incomplete network delineations[END_REF], the total network length that 1 Source: Conseil Général de l'Hérault governed the simulation process was not fixed. A thousand networks were simulated, and for each one, the targeted total network length was selected randomly and uniformly between a minimal value that corresponded to the total length of the invariant main downstream network, i.e., 10.8 km, and a maximal value that corresponded to the total cumulative length of the field units' boundaries, i.e. 220 km.

Hydrological modelling

The catchment hydrological behavior was simulated at the rainfall event scale using the physically based rainfall runoff MHYDAS model (Modélisation Hydrologique Distribuée des Agro-sytèmes-Distributed Hydrological Modeling of Agrosystems). This is a distributed model that considers the catchment as a series of interconnected geographical units by which infiltration/runoff partition and runoff routing are performed. Full details of the MHYDAS model description are available in [START_REF] Moussa | Hydrological modelling of a farmed mediterranean catchment[END_REF] and [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF]. The MHYDAS model has mainly been used in the context of farming [START_REF] Chahinian | Paramétrisation multi-critère et multi-échelle d un modèle hydrologique spatialisé de crue en milieu agricole[END_REF][START_REF] Charlier | Fonctionnement et modélisation hydrologique d'un petit bassin versant cultivé en milieu volcanique tropical[END_REF][START_REF] Gumiere | MHYDAS-Erosion: a distributed single-storm water erosion model for agricultural catchments[END_REF]Hallema et al., submitted;[START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF][START_REF] Tiemeyer | MHYDAS-DRAIN: A spatially distributed model for small, artificially drained lowland catchments[END_REF] and is available within the modeling platform Openfluid ® [START_REF] Fabre | OpenFluid: a software environment for modelling fluxes in landscapes[END_REF].

Landscape representation

We used the landscape representation of the MHYDAS hydrological model. The interconnected geographical units considered here were the field units (surface units) and the ditch network (linear unit), which was located on the boundaries of field units. The Geo-MHYDAS algorithm, running under GRASS GIS, was used to build a deterministic oriented tree topology between these irregularly shaped surface units and linear units, which allows the routing of simulated water flows across the landscape [START_REF] Lagacherie | Geo-MHYDAS: A landscape discretization tool for distributed hydrological modeling of cultivated areas[END_REF]. The Geo-MHYDAS algorithm uses as its inputs a DTM and the GIS layers of the ditch network and the field units. For each unit, the neighbor with the steepest slope is defined as its downstream unit.

In our case study, we used Geo-MHYDAS with the same 5 m DTM as the one used for network simulation. Geo-MHYDAS was thus used to define water flow paths of the surface units for each simulated network and for the actual network. The topological relationship between the field units and between the field units and the agricultural ditch drainage network were then computed for each simulated network and resulted in different topological configurations. In this process, water flow paths over the catchment (hillslopes and channels) were only modified by the ditch network spatial configuration.

Modeled hydrological processes

Over each surface (areal) unit, MHYDAS simulated the infiltration-runoff partition as a Hortonian process, while saturated runoff was neglected. The determination of the infiltration rate was based on equations from [START_REF] Green | Studies on soil physics part I: the flow of air and water through soils[END_REF] and [START_REF] Mein | Modelling infiltration during a steady rain[END_REF] that were adapted by Morel-Seytoux (1982) and allowed the calculation of the rainfall excess under variable rain conditions. It depended on the saturated hydraulic conductivity, which was a function of soil surface features that could vary in space and time. It also depended on the mean initial water content at the soil surface. The other parameters were the hydraulic properties of the hydrological units, such as the water content at saturation and the residual water content.

The rainfall excess function for each unit was converted to a surface runoff hydrograph by routing it to the proper outlet of the unit. The operation of converting excess rainfall into surface outflow was performed by a numerical convolution involving a unit hydrograph that was derived as a Hayami response function [START_REF] Hayami | On the propagation of flood waves[END_REF][START_REF] Moussa | Analytical Hayami Solution for the Diffusive Wave Flood Routing Problem with Lateral Inflow[END_REF]. This procedure was a function of the mean distances between the centers of gravity of adjacent hydrological units or between centers of gravity and the ditch network, the mean slope of the units and the roughness. If runoff was routed downstream via other surface units, it could be reinfiltrated to these surface units.

Once in the network that was considered as a set of linear reaches, water was routed to the catchment outlet using a diffusive wave equation. Over each reach, the parameters of the diffusive wave model were calculated using the length of the reach, its slope, the coefficient of roughness and the cross-sectional shape.

Groundwater was considered a compartment that received water that was infiltrated from surface units. Water could also be exchanged between the ditch network and groundwater according to the differences in water levels that were computed inside these compartments.

Parameterization strategy and simulation plan

Because our focus was the hydrological impact of the spatial configuration of the agricultural ditch drainage network, we considered a scenario in which the other sources of hydrological variability were deliberately removed. For this scenario, the corresponding parameters for the surface units (saturated hydraulic conductivity, water content at saturation, residual water content and roughness) and for the linear units (coefficient of roughness and cross-sectional size) were made invariant in the catchment. However, the values given to these parameters were selected so as to be coherent with the values that were calibrated in the previously mentioned studies of the same study area. They were also selected to simulate a coherent flood at the catchment outlet in comparison to the available outlet discharge data.

Rainfall was taken to be spatially homogeneous and was represented by a simple triangular rainfall of 50 mm in 4 hours to symbolize high and frequent rain events. Two other rainfall events, 30 mm in 2 hours and 60 mm in 4 hours, were used to test the sensitivity of the simulations to the rain. This set of three rain events represented events with a return period from one to two years in the area.

No interaction between surface and groundwater was introduced in the modeling because we focused on surface runoff and on the fast hydrological response of the catchment and interaction between the ditches and groundwater was a slower process. Moreover, during heavy rain events such as the one studied here, surface runoff was assumed to be the major component of runoff.

These simplified landscape and meteorological conditions allowed us to solely test the hydrological impact of the variability in the partitioning between diffuse and channelized flow paths and in the associated topographic parameters. Moreover, even if we did not aim to provide a validated case study, we used the model with a range of parameters and for a type of rain events in which the model has been calibrated and validated several times.

Evaluation of ditch network effects on hydrology

Classical indices

To compare simulated networks with the real network and to relate network runoff variability with network variability, we evaluated runoff indices (lag time, peak discharge and runoff volume) at three invariant locations within the catchment (points A-C, Figure 5). At the same three invariant locations, the total cumulative network length and the drained area were computed.

Need for new indices

To evaluate the variability in the runoff on the hillslopes (the overland flow), we defined an adequate geographical support at which the runoff could be computed. A 500 m × 500 m square grid was applied over the study area. Each cell of this grid was sufficiently large to include a piece of the landscape (individual fields were not) and could be also held invariant across simulations (subcatchments could not).

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

Inside a cell of this 500 m × 500 m grid, different hillslopes that were compounded of several field units were delineated during hydrological simulations. A hillslope was defined as a succession of field units whose overland flow exited the cell or entered the ditch network (Figure 4). The overland flow indicator of the cell was calculated as the mean of the maximum field unit (j) overland flow of all of the n hillslopes (i) inside the cell (Equations 1 and 2). For each cell, the cumulative network length was also computed.

Figure 4 Equation 1

Qmax hillslope i =max ( Qmax field j ) field j∈hillslope i j∈[1, number of fields inside hillslope i ] i∈[1, number of hillslopes inside the cell k ] Equation 2 Overland . flow . indicator cell k = 1 n × ∑ i=1 to n Qmax hillslope i hillslope i∈cell k k ∈[1, number of cells inside the catchment ]
3 Results

Agricultural ditch drainage network and water flow path variability

We simulated a thousand networks, which allowed us to explore the full range of potential cumulative network lengths (Figure 5). The invariant main downstream segments were conserved in both the actual network and simulated networks. Because of the simulation process, the mean length over the set of simulated networks was higher than the actual network, which was approximately half of the cumulative field boundary length. Some very different networks were obtained, including some that were very dense (Figure 5, IV) and others that were not very dense (Figure 5, II).

Figure 5 The modification of water flow paths by ditch networks implied modifications in subcatchment delineation and consequently the drained area (see the Roujan subcatchment

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422 defined at point B in Figure 5). To investigate this effect, we took as an example the evolution of the variability in the simulated drained area along the northern main invariant segment of the network (Figure 5, from point A to point C). At the start of the northern main invariant segment, the coefficient of variation of the simulated drained area was higher than 100 % (Figure 6). The coefficient of variation decreased rapidly with an increase in the mean simulated drained area. However, it was still high at point B (Roujan subcatchment) (20 %).

The coefficient of variation rapidly decreased to 5 % for a mean simulated drained area of 1.7

km². This high variability of the drained area emphasized the role of the ditch drainage network in modifying water flow paths. After the confluence with the southern main invariant segment, the coefficient of variation was low (decreasing to 2 %). The decreasing variability of the drained area along the main invariant segment could be explained by (i) the simulated drained area that became increasingly constrained by the catchment morphology and the invariant main segments near the outlet and (ii) the whole catchment boundary that was fixed for all of the simulations.

Discharge variability and scaling

We computed the hydrographs at three points along the northern main downstream segment of the network (Figure 7), including the starting point of the invariant main downstream network (Figure 5, point A), at 1000 m downstream (Figure 5, point B, corresponding to the Roujan subcatchment) and at the outlet, 5000 m downstream (Figure 5, point C, corresponding to the Bourdic catchment). The hydrographs exhibited a high variability in amplitude, but the shapes of the actual hydrographs were preserved across the simulated networks. For point A, the high variability that was observed was mainly related to the very high variability in the area that was drained at this point (Figure 6). In some cases, this point

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422 corresponded to a source, as no ditch was present upstream and it did not drain any field unit.

On the contrary, an upstream ditch network that drained several field units, could be branched upstream of this point for other simulations. For the two other points, the Roujan subcatchment and the Bourdic catchment, the variability in hydrographs remained high despite the smaller variability in the drained area. We hypothesized this property to be related to other network characteristics, especially network density.

Figure 7 We calculated lag times, peak discharges and total volumes from the 1,001 hydrographs. We focused on the sensitivity of these runoff indices to the variation in network lengths at two specific locations (Figure 8): at the outlet of a small upstream watershed, the Roujan subcatchment (the second hydrograph on Figure 6, point B on Figure 5) and at the outlet of the whole Bourdic catchment (the third hydrograph on Figure 6, point C on Figure 5). At the Bourdic catchment outlet, the range in peak discharge was -50 to +24 % compared with the value that was simulated for the actual network. The range in total volume was -55 to +14 %.

Furthermore, the peak discharge and total volume quickly increased with the network length and then reached a plateau. This plateau corresponded to a situation in which all of the fields units were connected to a ditch. Indeed, the denser the network, the greater the connectivity between the field units and the ditch network, which limited the re-infiltration to the fields.

The lag time exhibited a more complicated trend and less relative variability (-15 to +12 % in comparison with the value simulated for the actual network). First, there was a quick increase and then a slow decrease as the network length increased. We hypothesized that for short networks, only the runoff of the closest area to the outlet contributed to the flood and that the runoff of the more distant areas was reinfiltrated before reaching the network. For dense networks, distant areas contributed to the flood and thus the lag time was higher because of longer average flow paths. Then, when all areas contributed runoff and the network grew again, flow paths were increasingly channelized in the networks, which decreased the lag time because the mean celerity was higher in the network than in the field [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF].

Even if the catchment area was constant at the outlet, there was a variability in runoff for a given network length. This variability may have been caused by variability in the topographic parameters or the topology of the ditch networks.

The variability was more significant for the Roujan subcatchment scale. The lag time, peak discharge and total volume ranged from -18 to +90 %, -96 to +19 % and -71 to +29 %, respectively, in comparison with the simulated value for the actual network. An obvious reason was the fact that the drained area varied for the Roujan subcatchment outlet, but not at the Bourdic catchment outlet (Figure 6). However, there was also a clear trend for the peak discharge and total volume at the Roujan outlet: these variables increased as the network length increased. A trend was more difficult to identify for the lag time, but we observed a decrease with increasing network length.

Figure 8 Finally, we also investigated whether the relationship between the runoff metrics and the length of the ditch drainage network was dependent on the rain event. With two other triangular rain events, i.e., 30 mm in 2 hours and 60 mm in 4 hours, the type of relationship was not changed. An example of the peak discharge as a function of the ditch network length is shown in Figure 9.

Hillslope overland flow variability

Using the squared grid and the method described above (Figure 4), we observed how the overland flow indicator varied inside every cell grid on the catchment. Figure 10 shows an example of a cell that was located in the western part of the catchment, mainly in a zone with (approximately 2150 m in this cell) corresponded to a third of the maximum overland flow (6.9 l.s -1 against 21 l.s -1 ).For a given network length within a cell, a small amount of variability remained, which corresponded to the efficiency of the network to intercept overland flow. For example, if 2000 m of ditches were grouped into a small part of a cell or were parallel to the highest slope, they intercepted less overland flow than if they were well distributed throughout the cell and were perpendicular to the highest slope. We also noticed the position of the actual network, in that simulated networks with similar lengths produced more overland flow in the cell than the actual network.

As for the example cell in Figure 10 Figure 10 Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

Importance of the agricultural ditch drainage network on the hydrology of a small, cultivated catchment

Agricultural ditch drainage networks are known to influence the runoff of small, cultivated catchments [START_REF] Carluer | Assessment and modelling of the influence of man-made networks on the hydrology of a small watershed: implications for fast flow components, water quality and landscape management[END_REF][START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF]. However, the nature of this statement is dependent on the studied area and its actual network and is generally made through a comparison with a hypothetical natural network that is extracted from a DTM. In this study, the importance of ditch networks was confirmed in a more general context, even though it was still dependent on the catchment morphology and field geometry. Among the set of simulated networks, the variability in the area drained was very high, which reinforced the evidence that linear features alter overland flow-paths [START_REF] Duke | Incorporating ancillary data to refine anthropogenically modified overland flow paths[END_REF][START_REF] Gascuel-Odoux | Incorporating landscape features to obtain an object-oriented landscape drainage network representing the connectivity of surface flow pathways over rural catchments[END_REF]. Therefore, our study revealed how agricultural ditch drainage networks can modify delivery pathways, which is a major component of the hydrological connectivity of a catchment (Bracken & Croke, 2007).

The runoff variability was very high, both on the scale of the whole Bourdic catchment and the Roujan subcatchment. The peak discharge and total volume clearly increased with network density, but they also depended on other induced changes (e.g. topology and slope of the networks), which should be better quantified in the future. These results agreed with the study of [START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF], which concluded that agricultural ditch drainage networks accelerated runoff. However, we showed that this statement clearly depended on the drainage density. Concerning runoff on the hillslopes, we showed that the ditches that were present on hillslopes intercepted the diffuse water flow paths and rapidly decreased the overland flow.

However, above a certain network density, adding new ditches became less and less efficient for decreasing overland flow.

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

Our study also confirmed the predominant role of ditches in intercepting overland flow on the hillslopes in Mediterranean areas [START_REF] Martínez-Casasnovas | Soil erosion caused by extreme rainfall events: mapping and quantification in agricultural plots from very detailed digital elevation models[END_REF] and enabled the better quantification of the relationship between the density of the ditch network and the overland flow. This interception of overland flow thus explained the statistical effects of ditches in reducing erosion that was observed in the study by [START_REF] Paroissien | A regional-scale study of multidecennial erosion of vineyard fields using vine-stock unearthing-burying measurements[END_REF].

Finally, these results showed the effects of removing or adding ditches to a drainage network when no other field borders were considered when removing ditches. Conversely, in Verstraeten et al. (2002), replacing the ditches with grassed waterways was considered and was shown to reduce sediment delivery. Such an investigation was not considered here, and no obstacle to the runoff were considered between two given field units after removing a ditch.

Catchment extent effect

All of these findings were dependent on the extent of catchments. At any point of the catchment, a high variability in the network induced a high variability in runoff, but the variability in runoff was higher at the Roujan subcatchment scale than at the whole Bourdic catchment scale at any given network density. This fact was explained by the high variability in the area that was drained in the upstream part of the network, whereas the low variability in the area that was drained, the morphology of the catchment and the invariant main downstream network constrained the routing of the runoff and all served to limit the downstream variability (from approximately 2 km²).

Concerning the hillslope overland flow, the effect of scale was not investigated in this study.

However, the extent effect was considered to become significant with a decreasing cell size that would increase the overland flow indicator variability. A 500 m × 500 m square grid seemed to be appropriate because, in contrast to individual field units, it was sufficiently large

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

to include a piece of the landscape (which individual field units were not) and was sufficiently fine to enable the spatialization of overland flow in a small cultivated catchment with sufficient detail.

Hydrological efficiency

The actual ditch network seemed to be quite efficient in comparison with the simulated networks for draining floods (Figure 8) and limiting overland flow on the hillslopes (Figure 10). For instance, when looking at the peak discharge and the total volume of all of the points of the main invariant segments of the networks, the drainage that was provided by the actual network was almost always greater than the fitted value for the simulations, and the lag time was shorter (Figure 11). The actual network drained more water and accomplished that more rapidly. For the overland flow on the hillslopes, the overland flow indicator was almost always lower than that of the simulated networks; therefore, for a given length, the actual network intercepted more overland flow than the simulated networks. However, an optimum network should be not only more efficient than a network of similar length, but also a compromise between the network length and the hydrological efficiency because the maintenance of the ditch network requires a tremendous amount of work from farmers. To define such an optimum network, we must take into account the local conditions that influence the overland flow, such as land cover and soil properties. Therefore, we could not accomplish this task here. Finally, we question whether the actual agricultural ditch drainage network is optimal in all agricultural catchments or if this was the case in the Bourdic catchment only.

Figure 11 Current agricultural ditch drainage networks are the results of a sum of individual ditch digging efforts of hundreds of farmers over the last two thousand years [START_REF] Berger | Les fossés bordiers historiques et l'histoire agraire rhodanienne[END_REF],

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

and they continue to evolve. Therefore, we may wonder how this optimum can emerge from a sum of individual decisions. A survey of farmers would be useful to understand this finding and would perhaps highlight some guidelines for the creation of ditch drainage networks and cooperation between farmers, which could be further implemented in the simulation algorithm.

Limitations and uncertainties in results

This modeling study revealed the effects of modifications in the spatial configurations of ditch networks in an agricultural catchment, which can not be tested in the real world.

However, even if the model was previously used and validated in this area, we used this model in virtual conditions and outside its tested bounds (concerning drainage density).

Therefore, like in all studies testing scenarios, we must carefully interpret the results (Silberstein, 2006).

The results of this study should thus not be taken as universal, but the trends deserve consideration. In addition to these first limitations, the simplified hydrological case studied here does not allow for confidence in absolute values. Thus, this type of study should be extended to more complicated cases with various actual rain events and a realistically distributed land cover. Indeed, in our case study, each field unit generated the same amount of runoff, which would obviously be inaccurate in a real landscape. The consideration of patterns of runoff generation would probably modify the effect of the ditch networks. This effect would most likely be reinforced in areas with a high production of runoff and lessened in other areas. Ditch networks could thus modify the delimitation of active areas [START_REF] Ambroise | Variable 'active' versus 'contributing' areas or periods: a necessary distinction[END_REF], for example, by avoiding the runoff that is produced in a low-permeability area to be reinfiltrated to a permeable area that is located downslope.

Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

Another limitation concerns the parameterization of the ditch. First, neither infiltration into the ditch network nor ex-filtration were considered in this study, whereas groundwater recharge by channel infiltration may be important in Mediterranean and semi-arid areas (Dages et al. 2009;[START_REF] Ponce | Groundwater recharge by channel infiltration in El Barbon basin, Baja California, Mexico[END_REF]. If such an interaction was modeled, the network hydrology would probably have been much more complex, with patterns of infiltration and ex-filtration along the network, and the relation between network density and runoff may have been modified. Accordingly, we decided to focus on surface runoff in this study. Moreover, the calibration of the exchange function between a ditch network and groundwater is very delicate, especially with a variable network. Therefore, the interactions between surface and ground water in the ditches should be considered in future works, especially when studying rain events with lower intensities. The second limitation concerning the parameterization of the ditches involved their cross-sectional size and roughness which were considered invariant in our study. These parameters greatly affect surface runoff [START_REF] Nédélec | Réduction des transferts de crues dans les bassins versants agricoles fortement drainés : état des connaissances et des recherches[END_REF] and are known to be highly spatially variable, both in natural streams (O'Hare et al., 2010) and in small agricultural ditches [START_REF] Bouldin | Vegetative and structural characteristics of agricultural drainages in the Mississippi Delta landscapes[END_REF][START_REF] Crabit | A soft hydrological monitoring approach for comparing runoff on a network of small poorly gauged catchments[END_REF]. These limitations should therefore be addressed in future works, with the use of spatial models of these ditch parameters, as it exists for ditch cross-sectional sizes [START_REF] Bailly | Modelling Spatial Variability Along Drainage Networks with Geostatistics[END_REF].

To summarize, the results concerning the hydrological response of the network should thus be carefully interpreted. The results on overland flow in the fields seemed more robust because they were not affected by the simplification of the parameterization of the ditches.

Finally, the high variability that was observed here was predominantly a result of the high variability in the drainage density. However, some of this variability was improbable, especially the variability that resulted from including the very high drainage densities. submitted), the drainage density observed in the Bourdic catchment was one of the highest densities that has been observed. Therefore, the most interesting part of this study concerned low drainage densities for which hydrological sensitivity appeared to be the highest.

The coupling of a ditch network simulator with a hydrological model provided several interesting results and offered new prospects. Here, we used this method to investigate the impact of the spatial configuration of agricultural ditch drainage networks on the hydrology of a small catchment. We could imagine that this type of method could also be used to test certain scenarios of landscape arrangements or to optimize an existing network, as it has been done with land use patterns [START_REF] Newbold | A combined hydrologic simulation and landscape design model to prioritize sites for wetlands restoration[END_REF][START_REF] Seppelt | Optimization methodology for land use patterns using spatially explicit landscape models[END_REF]. For example, the minimization of the peak discharge at the outlet combined with the minimization of overland flow on hillslopes could be accomplished by branching or pruning the existing network. The optimization of the network could also concern reductions in erosion or pesticides impacts by coupling with appropriate models [START_REF] Gumiere | MHYDAS-Erosion: a distributed single-storm water erosion model for agricultural catchments[END_REF].

Conclusions

A ditch network simulator coupled with a distributed hydrological model allowed for a study of how runoff is related to the densities of an agricultural ditch drainage network in cultivated landscapes. New hydrological and network parameters were defined to deal with the variations in hillslope and subcatchment delineations that were induced by variations in the ditch drainage networks.

The importance of spatial configurations of agricultural ditch drainage networks in the alteration of water flow paths and in the control of runoff was highlighted, both in channelized flow-paths and on hillslopes. From our case study, we observed that variability of the spatial configuration and especially an increase of the ditch network density strongly influenced the subcatchment delineation and area, increased the drained volume and the peak discharge and decreased overland flow on hillslopes. However, these hydrological behaviors were only sensitive to network density for drainage networks with densities that did not exceed a certain threshold. This study also highlighted the efficiency of the actual network in comparison with the networks that were simulated by the ditch network simulator, which suggest an anthropogenic optimization process that should be further explored. Author-produced version of the article published in Hydrological Processes, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

Figure 2

 2 Figure2

Figure 6

 6 Figure 6

Figure 9

 9 Figure 9

  , a trend line could be fitted for each cell. To represent the exponential decrease of the overland flow indicator as a function of network length, a model in the form of overland.flow=a+b×exp(-c×network.length) was fitted. We used a non-linear least-squares algorithm for the estimation of the a, b and c parameters. The median coefficient of determination (R²) between the fitted and actual values of the overland flow indicator was equal to 0.96, which indicated a very good fit and an overland flow indicator that was closely dependent on the drainage density. The exponential decrease revealed that few ditches were efficient in greatly limiting overland flow and that more and more ditches were needed to further reduce the overland flow.

Figure 1 :

 1 Figure 1: Linear features in cultivated landscapes. Ditches and roads channel water, whereas embankments act as barriers to water flow.

Figure 5 :Figure 7 :

 57 Figure 5: Variability of the spatial configuration and length of the drainage networks. Both small and very dense networks were simulated. We also observed variability in the delineation of subcatchments.

Figure 8 :

 8 Figure 8: Relationship between runoff metrics and network lengths. The runoff metrics were correlated with the drainage network length at both the subcatchment and at the catchment scales. The variabilities in network metrics were higher at the subcatchment scale.

Figure 9 :Figure 11 :

 911 Figure 9: Relationship between peak discharge and network length for three different rain events.

  Author-produced version of the article published in HydrologicalProcesses, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422 slopes greater than 10 %. The overland flow indicator decreased exponentially as the cumulative ditch network length increased and trended toward an asymptote across the highest drainage densities. Approximately 25 % of the maximum drainage density

  Author-produced version of the article published in HydrologicalProcesses, 2012, 26(22), 3393-3404. The original publication is available at http://onlinelibrary.wiley.com DOI:10.1002/hyp.8422

However, their inclusion provided the advantage of defining what was possible. Compared with the large survey that was realized in the Hérault French département

(Levavasseur et al.,

Acknowledgement

This work was partly supported by the French Languedoc-Roussillon région, via a Ph.D. thesis grant. The authors also thank Adeline Bellet, an M.Sc. student, who helped with field surveys.