Minh-Ngoc Vu

Ahmad Pouya

Darius M Seyedi

Modelling of steady-state fluid flow in 3D fractured isotropic porous media: Application to effective permeability calculation

Keywords: Fractured porous media, steady-state fluid flow, effective permeability, Upscaling, Integral Equation, Poiseuille's flow

come     Modelling of steady-state fluid flow in 3D fractured isotropic porous media: Application to effective permeability calculation

INTRODUCTION

In the last two decades, there has been an intensive investigation in the fluid flow analysis in cracked porous media or fractured rock masses, due to the modelling requirements for important industrial applications such as: underground radioactive waste repositories, natural oil/gas recovery, geological CO 2 storage, geothermal energy, etc. However, determining effective permeability for fractured rock formations is one of the most challenging problems that are faced by geomechanical researchers and geotechnical engineers. Effective continuum is a simple and efficient approach to study the macroscopic properties of fractured and heterogeneous media by taking over a representative elementary volume. Discrete Fracture Network (DFN) model has been a powerful approach to simulate the flow in the natural fractured reservoirs for a long time [START_REF] Long | A model for steady fluid flow in random threedimensional networks of disk-shaped fractures[END_REF][START_REF] Cacas | Modeling fracture flow with a stochastic discrete fracture network: calibration and validation[END_REF][START_REF] Pouya | Permeability of 3D discontinuity networks. New tensors from boundary-conditioned homogenisation[END_REF][START_REF] Jafari | Effective Fracture Network Permeability of Geothermal Reservoirs[END_REF]. DFN is introduced in several commercial softwares as UDEC (Itasca) [START_REF]Universal Distinct Element Code (UDEC)[END_REF], Fraca [START_REF]Fractured Reservoir Characterization and Modelling (FracaFlow)[END_REF], etc. Nevertheless, the matrix is supposed to have negligible porosity and permeability in DFN models. A common approach to overcome this limit is dual porosity/dual permeability model that was first proposed by Barenblatt et al. [START_REF] Barenblatt | Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[END_REF]. In this framework, fracture and porous matrix are considered as two interactive continua. This approach has been widely used for modelling fluid flow in fractured materials. A brief review of several flow models can be consulted in Jafari and Babadagli [START_REF] Jafari | Effective Fracture Network Permeability of Geothermal Reservoirs[END_REF]. Different numerical techniques, such as finite element method (FEM) [START_REF] Baca | Modeling fluid flow in fractured porous rock masses by finite element techniques[END_REF][START_REF] Ghabezloo | Numerical upscaling of the permeability of a randomly cracked porous medium[END_REF][START_REF] Woodbury | Lanczos method for the solution of groundwater flow in discretely fractured porous media[END_REF][START_REF] Lesinigo | A multiscale Darcy-Brinkman model for fluid flow in fractured porous media[END_REF] and finite volume method (FVM) [START_REF] Castaing | Taking into account the complexity of natural fracture systems in reservoir single-phase flow modelling[END_REF][START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF][START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF][START_REF] Grillo | Simulation of density-driven flow in fractured porous media[END_REF][START_REF] Mourzenko | Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations[END_REF] have been used to resolve fluid flow equations in fractured porous media. Modelling of the steady-state flow in a fractured porous media under a far-field condition constitutes a key problem for effective permeability determination. In this work, we investigate numerically the flow around several fractures embedded in an infinite body described by a uniform gradient pressure at infinity. The fracture is modelled by a zero thickness surface with no pressure jump between two faces; i.e., having an infinite transverse conductivity. In order to upscale effective permeability, it is interesting to consider the solution of flow around an isolated fracture. This problem can be derived by using the Darcy's ellipsoidal inclusion model [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF][START_REF] Dormieux | Approche micromécanique du couplage perméabilité endommagement[END_REF][START_REF] Barthélémy | Effective Permeability of Media with a Dense Network of Long and Micro Fractures[END_REF]. Nonetheless, Poiseuille's law is a more common model to describe the flow in a single fracture in the physical rock area. A comparison between these two models can be seen in Pouya and Vu [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF]. Potential solution based on singular integral equation is proposed recently to study the steady-state flow in an infinite fractured porous medium with homogeneous permeability matrix. The 2D problem is studied by Liolios and Exadaktylos [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF]; Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]; Pouya and Vu [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF]. In these studies, the effective permeability is computed analytically by Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] and numerically by Pouya and Vu [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF]. The potential solution was extended to 3D problem by Pouya [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF], which reveals the pressure field in porous matrix as a function of the infiltration through the curved fracture surfaces. Therefore, this advance of mathematical formulation allows us to decrease dimension of the discretised domain: 3D (pressure field in the whole body) to 2D (infiltration through the fracture surfaces). Pouya and Vu [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF] used this potential solution to establish a semi-analytical solution for a superconductive elliptical discshaped fracture governed by Poiseuille's law flow. Moreover, a closed-form solution for a superconductive circular disc fracture is derived in this work [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF]. The obtained solution is then used to estimate the effective permeability. However, any analytical solution for the general case of the fracture with a finite conductivity does not exit. Hence, numerical resolution is naturally imposed. This paper focuses on singular integral equation method to resolve the potential equation in 3D fractured porous media. Some examples with several intersecting fractures are presented. Numerical solution for a single fracture obtained in this way permits us to establish a semi-analytical description of infiltration through the fracture as a function of fracture conductivity, matrix permeability and fracture geometry. Then, Mori-Tanaka and self-consistent schemes are employed to estimate the effective permeability of a fractured porous media by using this semi-analytical solution. The superconductive fracture is always used in this work as a reference problem to make reader confident on the results.

Notations

In the following, light-face (Greek or Latin) letters denote scalars; underlined letters denote vectors, bold-face letters designate second rank tensors or double-index matrices. The scalar product of two vectors a and b is labelled as a.b. For second rank tensors, the tensor transposed from A is denoted A T , the matrix product is labelled as AB and the determinant as |A|. The operation of A on a is labelled as A.a. The convention of summation on repeated indices is not used for Latin indices (i, j, k…) that number surfaces, lines, etc., and are noted indifferently as subscript or superscript. ∇ represents the gradient and ∆ the Laplace operator for a scalar field and (∇.) the divergence for a vector field.

GOVERNING EQUATIONS

Considering a reference volume Ω of an infinite three-dimensional porous media having isotropic permeability matrix k=kδ and containing a family of fractures numbered by j and noted Γ j (Fig. 1). Each fracture surface Γ j is defined by a smooth function z j of the curvilinear parameters {s α } (α=1,2). Two or more fractures can intersect along an intersection line and three or more fractures can cross at an intersection point. The intersection lines and the border of fracture surfaces ∂Γ constitute a set of singular points denoted by S [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF]. The set of points on the fracture surface is denoted by Γ. Therefore S ⊂ Γ and rock matrix corresponds to Ω-Γ. A far-field linear pressure defined as p ∞ (x)=A.x is applied to the domain Ω, with A a constant pressure gradient.

The fluid flow in the matrix is governed by Darcy's law:

∀x∈Ω-Γ ( ) ( ). ( ) v x x p x = - ∇ k (1)
where, v(x) and p(x) are local fluid velocity and pressure, respectively.

∀x∈Ω-Γ . ( ) 0 v x ∇ = (2)
As mentioned by Guéguen and Palciauskas [START_REF] Guéguen | Introduction to Physics of Rocks[END_REF], Poiseuille's law is the most common model for flow in a single fracture. This law establishes a linear relation between the flux vector and the in-plane pressure gradient in the fracture as follows:

∀s∈Γ ( ) ( ). ( ) s q s s p s = - ∇ c (3)
where, c is the hydraulic conductivity and consists of a tensor of 2x2 dimension in the plane of fracture and is assumed to be isotropic (i.e., c = cδ). By neglecting roughness effect, the hydraulic conductivity of a fracture c is deduced by the classical cubic law.

Hence c depends on the fracture aperture e and the fluid dynamic viscosity µ, as c=e 3 /(12µ). However, real rock fractures have rough walls, variable aperture and asperity areas. To simplify the problem, a single average value, called hydraulic aperture e * , is usually employed to characterize the mechanical aperture e of an individual fracture. The empirical relations between these two variables can be found in a great amount of works in the literature. A review of these models is available in Zhou et al. [START_REF] Zhou | Flow-stress coupled permeability tensor for fractured rock masses[END_REF].

The in-plane pressure gradient is defined by the flowing formulation:

1 2 1 2 ( ) s p p p s t t s s ∂ ∂ ∇ = + ∂ ∂ (4)
and, (t 1 ,t 2 ) defined by t α =∂z/∂s α are two perpendicular unit vectors tangent to Γ (Fig. 2)

n z t 1 t 2 s 2 s 1 D Γ Fig. 2
Curvilinear parameters [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF] The fracture-matrix mass exchange at a regular point on the fracture Γ-S (excluding the intersection points) is widely described in the literature [START_REF] Barenblatt | Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[END_REF][START_REF] Adler | Porous Media: Geometry and Transports[END_REF][START_REF] Aifantis | On Barrenblatt's problem[END_REF][START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF][START_REF] Mourzenko | Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations[END_REF]. This law generally has been obtained by considering mass balance in a small volume D around a point z. Using the scheme presented in Fig. 3a, the mass conservation reads [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF]:

( ) . ( ) . ( ) 0 s v z n s q s + ∇ = (5)
where, n is the normal unit vector to the fracture surface, pointing from Γ -to Γ + . v -and v + are the fluid velocity in the matrix at the point z of lower and upper fracture surface, respectively. ( )

( ) ( ) v z v z v z + - = -
designates the discontinuity or jump of fluid velocity across the fracture. ∇ s .q(s) represents the surface divergence. If the curvature effect of surface is neglected [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF], the surface divergence reads:

1 2 1 2 ( ) s q q q s s s ∂ ∂ ∇ = + ∂ ∂ (6)
Mass balance at intersection lines between several fractures surfaces is an important issue in studying the fluid flow in the fractured porous media. Pouya [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF] presented a rigorous way to obtain the closed-form formulation of mass balance at a singular point. This relation is similar to one that can be found for a fracture network in an impermeable medium.

. 0

j j j q m = ∑ (7)
where, q j is the infiltration vector on Γ j and m j tangent normal vector of the intersection line on Γ j pointing outward from the domain as seen in Fig. 3b. The pressure field p(x) in a fractured porous domain, under a far-field constant pressure gradient, must satisfy the constitutive equations (1), (3) and mass balances equations (2), ( 5), [START_REF] Barenblatt | Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[END_REF] and boundary conditions. For fractures with infinite transversal transmissivity (no pressure jump between their two sides), Pouya [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF] proposed a general pressure potential solution as a function of the infiltration through the fractures such as:

q(s) v - D + D - n z v + m D ∂D (a) S + S - S j - γ 1 S j - + γ P L γ L P j γ S j + γ 1 P j+ γ V j Γ 1 S j + + γ S j - γ L γ L m m j m m V j Γ P j γ S j + γ j q (b)
( ) 3 ( ). ( ) 1 ( ) ( ) 4 ( ) j j j j j q s x z s p x p x ds k x z s ∞ Γ - = + - ∑ ∫ π (8)
For the effective permeability estimation purposes, the numerical calculation is only applied for plane surface factures in the present work. Nevertheless, equation [START_REF] Baca | Modeling fluid flow in fractured porous rock masses by finite element techniques[END_REF] is valid for the curved intersecting fracture surfaces with arbitrary forms.

Introducing equation ( 3) into (8), a hypersingular integral equations is obtained, in which the pressure field on fractures is the main unknown. The proposed flow model makes an important advantage by reducing a three-dimensional problem to two-dimensional one. Hence, 2D mesh is only generated on the fracture surfaces in the first stage to obtain the pressure and the infiltration solutions on the fractures, and then, these results are used again to compute the pressure field in the whole domain in the second stage. Unlike FEM and FVM, which pass by volume discretization and require 3D meshes in the whole space.

NUMERICAL RESOLUTION

FEM is very efficient for solving finite-domain problems, but it is less convenient to deal with infinite domain problems. Integral equation method or boundary element method is an efficient alternative that has been extensively used for solving infinite or half-plane problems [START_REF] Brebbia | Boundary Elements: An Introduction Course[END_REF][START_REF] Brebbia | Boundary Element Techniques Theory and Applications in Engineering[END_REF]. This section describes a numerical procedure to solve equation ( 8) based on collocation method. The principle resolution of integral equation by collocation method can be found in Bonnet [START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF]. A robust methodology is proposed in this work to achieve fast numerical calculations. The numerical procedure is presented for solving a general problem to model the flow in an infinite isotropic medium containing several intersecting curved fracture surfaces under a uniform far-field pressure gradient. Some results of pressure field on the fractures are illustrated for many intersecting fractures or a curved surface at the end of this section to demonstrate the capacity of the developed numerical code.

Geometry discretization and variables interpolation

A mesh is generated on the surface of all fractures. It is constituted by N nodes and M triangular linear elements numbered m and denoted E m . This mesh must be consistent at the intersection lines between the fractures. By replacing the infiltration by equation ( 3), an approximated form for equation ( 8) can be obtained:

1 ( ) ( ) ( ) 4 M m m c p x p x I x k π ∞ = ≈ - ∑ (9)
where, I m is an elementary integral:

( ) ( ) 3 . ( ) ( ) ( ) m m E p z x z s I x ds x z s ∇ - = - ∫ (10)
The choice of the interpolation function for the elements at the boundary of fractures is based on the velocity profile in this area. Regarding one-dimensional singular integral equation over a finite interval, many authors [START_REF] Muskhelishvili | Singular Integral Equations[END_REF][START_REF] Erdogan | Approximate Solutions of Systems of Singular Integral Equations[END_REF][START_REF] Martin | End-point behaviour of solutions to hypersingular integral Equations[END_REF] characterized the solution profile at the extremity points, which varies as r 1/2 in vicinity of these points. Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] derived the problem of fluid flow around an isolated superconductive crack line located at the interval [-1,1] on the Ox-axis and under a constant gradient of pressure at infinity. A closed-form solution of flux is found, which varies as 2 1 r in the fracture. Pouya and Vu [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF] used the interpolation of discharge q(s) corresponding to a variation according to s 1/2 for the general case of Poiseuille's fracture, with s the distance from extremity point. In two dimensional, Wolfe [START_REF] Wolfe | Eigen/functions of the integral Equation for the potential of the charged disk[END_REF] and Boersma and Danicki [START_REF] Boersma | On the solution of an integral Equation arising in potential problems for circular and elliptic disks[END_REF] considered the problem of the electrostatic potential due to a charged elliptic disk. The solution also varies as r 1/2 closed to disk border, with r the distance from the disk border. Recently, Pouya and Vu [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF] analytically derived the problem presented in this work for a single superconductive circular disc fracture. They showed the infiltration solution varies as As a result, three types of element are distinguished: Type 0, 1 and 2 correspond to the element that has no node, one and two nodes on the border of fractures respectively (Fig. 4). Element of type 0 is considered as isoparametric, i.e. a linear interpolation of pressure is used that leads to a constant flux q on that elements. On the other hand, inspiring from closed-form solution of two-dimensional and three-dimensional problems, the interpolation function of flux vector is chosen for element of types 1 and 2, which varies as r 1/2 in the perpendicular direction of the facture boundary and be constant in the other direction in each element of those types as depicted in Fig. 5. Interpolation functions of pressure on the elements of type 1 and 2 can be deduced from equation ( 3) and its formulations in global coordinates systems is presented in Appendix A.

To simplify the notation, ( ) k x ψ is used to show the interpolation functions of three types of elements throughout this text.

Construction and resolution of linear system of Equations

The collocation method consists in forcing the integro-differential equation ( 9) to be verified exactly at collocation points x c . M collocation points are chosen as the set of barycenters of M triangular elements and are numbered as the elements. Therefore, the discretization of equations, for collocation point i c

x on the element i, is performed by transferring the variables interpolations:

1 ( ) ( ) ( ) 4 M i i i c c c m m c p x p x I x k π ∞ = ≈ - ∑ (11)
such as:

( ) ( )

3 3 3 1 1 1 ( ). . ( ) ( ). . 4 ( ) m i m M c k k i i i c c k k i k m k E c z p x z s c x p A x ds k x z s = = = ∇ ψ - ψ ≈ -π - ∑ ∑∑ ∫ (12)
This equation is first written for M collocation points, and then, assembly operation based on two nested loops over all collocation points and all elements lead to:

. = H P G [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF] where, G(Mx1) is the column of pressure prescribed at infinity p ∞ (x) at M collocation points, P(Nx1) the column of nodal pressures or unknowns, and H(MxN) the matrix of coefficients obtained by elementary integration calculation explained in the next section.

It should be noted that the system of Equations ( 13) consists of M equations and N unknowns with M>N. Indeed, most of meshes built on the surface of a set of fractures contain more elements than nodes. This situation is always seen for fine meshes. However, the approximated solution could be computed by a least squares fitting method minimizing |H.P -G| 2 .

Elementary integral calculation

Equation ( 12) contains an elementary integral: (

m i m c k k mk i E c z p x z s I ds x z s ∇ ψ - = - ∫ and 3 1 m mk k I I = = ∑ ) ( ) 3 ( ). . ( ) ( )
This integral is singular if

i c m
x E ∈ .i.e. i=m and regular if

i c m
x E ∉ or i≠m. The singular elementary integrals are evaluated in Cauchy Principal Value meaning [START_REF] Guiggiani | The evaluation of Cauchy principal value integrals in the boundary element method-A review[END_REF]. In the boundary element method field, computation of regular elementary integral is usually based, as well as in finite element method, on the use of Gaussian quadrature formulas. However, ordinary numerical method would not yield reliable and consistent solutions for evaluating singular elementary integrals since the singularities arise from the kernels of integrals. Several powerful schemes are applied to compute accurately these integrals such as: integration by parts [START_REF] Bonnet | Methodes des Equations Integrates Regularises en Elastodynamique[END_REF], Limit process [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF], and transformation method [START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF][START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF][START_REF] Karami | An efficient method to evaluate hypersingular and supersingular integrals in boundary integral quations analysis[END_REF]. The numerical difficulty appears when the collocation point is too close to the element, so the integral is not strictly singular but nearly singular. Boundary element method researchers have proposed various numerical schemes or semi-analytical approaches to derive such nearly singular integrals. A brief review of some techniques can be found in Liu et al. [START_REF] Liu | Nearly singular and hypersingular integrals in the boundary element method[END_REF].

In this work, the elementary integral may be computed by different ways, namely: analytical, numerical and semi-analytical depending on the relative position of collocation point x c regarding the element E and the type of element (i.e. interpolation function). .

Analytical scheme

Analytical procedure is developed to calculate the singular elementary integral over element of type 0 (Fig. 6a). The treatment technique of singularity is based on the introduction of a polar coordinate system centred at x c in element E:

c z x = + ρφ (15)
where, (cos ,sin ) φ = θ θ . Linear interpolation of pressure over element leads that the gradient of interpolation functions is constant, thus, we suppose ( )

( ) k k z const ∇ ψ = ζ = .
Element E is divided into four parts: a small circle C of radius ε centred at x c and three sub-triangles T l (l=1,2,3).

Rewriting integral [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF] in polar coordinate yields:

3 3 1 1 l m mk k k k l C T d d I p d d = =   ρ ρ = ζ φ θ + φ θ     ρ ρ   ∑ ∑ ∫ ∫ (16) then, 2 1
( ) 2 [START_REF] Pouya | Permeability of 3D discontinuity networks. New tensors from boundary-conditioned homogenisation[END_REF] 3

1 1 0 0 l l l m mk k k k l d d I p d d θ ρ θ ε π = = θ ε   ρ ρ = ζ φ θ + φ θ     ρ ρ   ∑ ∑ ∫ ∫ ∫ ∫ (17)
where, the parameter of each sub-triangle T l (l=1,2,3) in equation ( 17) can be found in Fig. 7. As expected, the first integral in this Equation is a singular integral which is vanishing due to

∫

. The second integral is a regular one and can be integrated explicitly. Indeed, ε appears in the result of integral over each sub-triangle and vanishes after summing three integrals.

Let us consider ϒ as the plane that contains element E. It is interesting to notice that this analytical scheme is also valid for the case in which the collocation point x c is on plane ϒ and is not in the element E of type 0 (Fig. 6b). This procedure provides a fast numerical calculation for evaluating the effective permeability of fractured porous media as demonstrated in section 0.

Semi-analytical scheme

Singular elementary integrals cannot be computed explicitly over the element of types 1 and 2 due to complicated nature of interpolation function in equations ( 49), (50). A semi-analytical procedure is developed to evaluate these kinds of integrals, which is based on the transformation to polar coordinate system centred at collocation point x c such as one presented in the aforementioned scheme. In view of singularity treatment, the integration of variable ρ is performed analytically. In this manipulation, the singularity of ρ is disappeared. And then, the integration of variable θ is regular and can be integrated by Gaussian quadrature in one-dimensional.

In addition, the nearly singular elementary integral is also computed by this semianalytical procedure in which the polar coordinate system is centred at the projection x pc of collocation point x c on element E (Fig. 6c). Supposing h as the distance between x c and element E, integral [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF] for element of type 0 in polar coordinate reads: As analysed in the analytical scheme in equation ( 17), only the second integral in this formula contributes to final result in which the integration of ρ can be manipulated explicitly and then a numerical integration is employed for variable θ. Unlike the general numerical procedure of boundary element method, the nearly singular integral is evaluated by some transformation technique for concentrating the Gauss's points in the vicinity of x pc that costs many computational time and memory of PC.

( ) ( ) 2 1 ( ) 2 2 2 3 3 3/2 3/2 2 2 2 2 1 1 0 0 l l l m mk k k k l d d I p d d h h θ ρ θ ε π = = θ ε   ρ ρ ρ ρ   = ζ φ θ + φ θ   ρ + ρ +   ∑ ∑ ∫ ∫ ∫ ∫ (
The importance of correct evaluation of singular and nearly singular elementary integral takes full advantage for the good conditioning of linear system. Especially, the application of this numerical method in fractured porous rock formation in which the ratio between fracture conductivity and permeability of rock is normally considerable.

Numerical scheme

Let us consider x cp as the projection of x c on ϒ. Gaussian quadrature rule is an efficient method for computing the elementary integral in the case x c is not on the plane ϒ and x cp is not in E (Fig. 6d). Referring to the theory of numerical integration [START_REF] Davis | Methods of numerical integration[END_REF], the approximate value of integral ( 14) reads:

( ) ( )

3 3 1 1 ( ) . i Ng Gj c Gj k m m k j i k j c Gj x x x I p w x x = = ∇ ψ - = - ∑ ∑ (19)
where, Ng is the number of Gaussian point, x Gj the Gaussian point and w j weight of number j. This scheme can be also used to evaluate regular elementary integrals over the elements of type 1 and 2 where the collocation point x c is on the plane ϒ

VALIDATION OF NUMERICAL SOLUTION

The limit case of a superconductive fracture is used as a benchmark problem for the numerical method in this work due to the fact that the closed-form solutions are well known for this reference case. In this view, two test-cases with different geometries are proposed.

Single superconductive disc fracture

The first case, an isolated superconductive fracture (c→∞) is modelled by a zerothickness disc of radius R in an infinite porous matrix under the limit condition p ∞ (x)=x. This disc is placed on the plane XY as depicted in Fig. 8. Pouya and Vu [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF] presented a rigorous procedure to derive the explicit solution for a superconductive fracture (i.e. the pressure is constant in the fracture surface). The analytical solution of flux in the fracture disc is expressed in equation ( 20) and its contour is shown in Fig. 9.

Let us now apply the numerical calculation for a single superconductive fracture (c>>1 for numerical modelling). The numerical solution of pressure obtained by singular integral equation is constant in superconductive fracture. Regarding the flux solution, a perfect agreement between analytical solution and numerical solution is also observed (Fig. 10). A relative error is then calculated to compare numerical and analytical results [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF], where, q n , q a are infiltration solutions obtained by numerical and analytical ways, respectively. The maximum value for the error is E rmax (q x ) = 0,06%.

( )

3 3 2 3 2 2 2 2 2 3 2 2 2 8 ˆ( ) 3 2 2 3 k q Rr R R r r R r R R r A r -      ξ = -+ - + -+ - - ξ ⊗ ξ       π       T δ δ δ δ (20)
where, . q(r,θ)

r = ξ ξ , 1 ˆr ξ = ξ , T = δ δ δ δ -e z ⊗ e z .
-q ax (r,0) -q nx (r,0) -q ax (r,π/4) -q nx (r,π/4) -q ax (r,π/2) -q nx (r,π/2) q ay (r,π/4) q ny (r,π/4) Fig. 10 Analytical vs. numerical solutions of infiltration q(r) for different value of θ. 0 ≤ r ≤ R and q ax , q ay , q nx , q ny are analytical and numerical infiltration field on the fracture in x and y directions

Single superconductive band fracture

In the second case, a single superconductive fracture (c→∞) is modelled by a zerothickness band of finite width L and of infinite length Λ (numerically Λ>>L) in a 3D infinite porous medium of isotropic permeability k=kδ (Fig. 11 left). Pressure field p ∞ (x)=y is imposed at infinity. Theoretically, the fluid flow in every plane being perpendicular to x-axis is similar. As a consequence, this 3D configuration is equivalent to the fluid flow in an infinite plane (plane YZ) containing an isolated fracture line of length L under the limit condition p ∞ (x)=y (as seen in Fig. 11 right). For this 2D configuration, a closed-form solution of flux in fracture line is given by Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]. From a numerical point of view, a band fracture with Λ/L=5 is modelled. The numerical solution of pressure field in fracture allows computing the flux in y-direction that passes at the middle of fracture. Good agreement between the current approach and the benchmark solutions is observed in Fig. 12. These two examples of validation show that the developed numerical procedure provides results with high fidelity and accuracy.

APPLICATION TO FLOW THROUGH THREE-DIMENSIONAL FRACTURED POROUS MEDIA

In order to demonstrate the capacity of the developed numerical procedure to simulate the flow through fracture networks, two application examples are presented in this section. The first one consists in modelling fluid flow through several intersecting circular fractures. In the second example, a curved fracture surface with a constant conductivity embedded in an infinite matrix is considered. In the first case, a consistent mesh is generated on the circular fractures as observed in Fig. 13. Fig. 14 presents the obtained pressure field in the fracture network. As observed, the pressure seems be continuous on the intersection lines between fractures. The masse balance ( 7) at these lines can be verified and is well satisfied as it was also checked in our previous work on 2D problem [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF]. In the second one, a curved fracture surface is modelled as a piece of a spherical surface (Fig. 15a). The numerical solution of pressure on this fracture is illustrated in Fig. 15b. However, we do not intend to analyse these results since it is not the main objective of this paper and we hope return to this subject in an application work.

APPLICATION TO EFFECTIVE PERMEABILITY ESTIMATION

Contrary to the case of superconductive disc fracture, the case of disc fracture having a finite conductivity has not an analytical solution. Numerical computation based on singular integral method can deal with this problem. The numerical procedure described in section 3 is applied for a single circular fracture with different conductivity values.

Obtained numerical solutions allow us to determine a semi-analytical solution of infiltration through this fracture and effective permeability of fractured porous media as a function of fracture conductivity, matrix permeability and geometry variables.

Semi-analytical solution of flux on fracture

An infinite body with isotropic permeability k=kδ containing a single fracture disc of radius R and of constant conductivity c is considered (Fig. 8). The body is submitted at infinity to a constant gradient A. Equation ( 8) is written for this problem such as:

( )

( ). ( )

. 4

p z x z c p x A x dz k x z Γ ∇ - = -π - ∫ (22)
Let us introduce following dimensionless variables:

x X R = , 4 c kR λ = π (23)
And the notation, ( )

( ) p x P X R A λ = , ( )
( ) p x A P X λ ∇ = ∇ (24)
By supposing Γ* as a circle of unit radius centred at the origin of coordinates system, equation

. P X X Z A P X X dZ A X Z λ λ Γ ∇ - = -λ - ∫ (22) reads: ( ) ( ) ( ) * 3 .
Numerical procedure described in section 3 is employed to resolve equation ( 25) for calculating the pressure field P λ (X) corresponding to different values of λ. The numerical result for λ = 2 is shown in Fig. 16. Evaluation of integral flux field over the fracture surface gives:

* 2 ( ) ( ) ( ) Q x c p x dx cR A P X dX λ Γ Γ = -∇ = - ∇ ∫ ∫ (26)
By defining, The value obtained numerically for P λ (X) allows describing with a great precision the relation between two variables Q * (X) and λ by the following formulae as plotted in Fig. 17.

* * 2 2 ( ). ( ) . ( ) Q x A Q X c A P X dX R A λ Γ = = - ∇ ∫ (27)
* 3 1 4 3 c Q π = λ + (28)
Despite the fact that the above expression is obtained numerically by curve fitting, it reveals to be extremely accurate. Equation [START_REF] Aifantis | On Barrenblatt's problem[END_REF] provides a relationship between integral of flux field over the circular disc fracture ( ) Q x and the problem parameters, namely: fracture conductivity c, matrix permeability k, and fracture geometry R:

3 16 ( ) 16 3 3 kcR Q x A c kR = - + (29)
This equation can be verified by its asymptotic behaviour when λ→∞ such:

3 16 lim ( ) 3 c kR Q x A →∞ = -
, which is according to close-formed solution for superconductive circular disc fracture given by Pouya and Vu [20]. Furthermore, this result is consistent with one obtained by Shafiro and Kachanov [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF] for the superconductive flattened ellipsoidal inclusion with Darcy's flow. This property reveals a high precision of the semi-analytical solution [START_REF] Brebbia | Boundary Elements: An Introduction Course[END_REF]. Therefore, this result can be used confidently in the following to calculate the effective permeability of a fractured porous media.

In order to determine the equivalent permeability tensor of fractured porous media, equation ( 29) could be expressed for the general case of an arbitrary circular fracture in space such as: where, T is the projection operator on the plane of fracture given by the following expression in which n represents the unit normal of the fracture plane:

= n n -⊗ T δ δ δ δ (31)

Upscaling of effective permeability

Let us consider a set of fractures Γ m embedded in a reference homogeneous volume Ω of an infinite 3D porous medium having isotropic permeability k = kδ. By considering the fracture as the limit case of a thin permeable layer, the average velocity in Ω reads [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF][START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF][START_REF] Mourzenko | Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations[END_REF][START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF]:

1 m m V vd qd Ω Γ   ≡ Ω + Γ   Ω     ∑ ∫ ∫ (32)
A linear pressure p(x) = A.x is applied at the boundary ∂Ω. Using the divergence theorem for the pressure function in a domain containing des discontinuities with the fact that p =0 on the fractures, allows the determination of the average pressure gradient as follows:

1 G pd A Ω ≡ ∇ Ω = Ω ∫ (33)
Introducing (1) into the first term of V in (32) yields:

1 . vd A Ω Ω = - Ω ∫ k (34)
In addition, the effective permeability k is yielded by V = -k .G, hence, ˆf = k k + k [START_REF] Wolfe | Eigen/functions of the integral Equation for the potential of the charged disk[END_REF] where, k f is the contribution of fractures to the effective permeability defined as:

1 .

m f m qd A Γ Γ = - Ω ∑ ∫ k (36)
According to the above-mentioned analyses, the semi-analytical solution expressed in equation ( 30) could be used to calculate the effective permeability in equation [START_REF] Wolfe | Eigen/functions of the integral Equation for the potential of the charged disk[END_REF] when the domain Ω tends to infinity. For a weak fracture density, dilute Mori-Tanaka scheme for estimation of k can be used. In this scheme, the interaction between fractures is neglected and each fracture is subject to the same far-field pressure, unperturbed by the presence of other fractures. In this approximation, the effective permeability tensor k is expressed as the sum of the matrix permeability and additional permeability due to each fracture k f :

ˆf m m k = + ∑ k k δ δ δ δ (37)
Assuming that each family of fractures Γ m is characterized by a mean radius R m , an orientation n m (unit normal of fracture), a conductivity c m and a density (number of fracture center per unit volume) ρ m , the expression of f m k will be obtained as follows:

( )

3 16 16 3 3 m m m m f m m m m kc R n n A c kR = ρ -⊗ + . k δ δ δ δ (38)
Equations ( 37)-( 38) cover arbitrary orientations of fractures and constitute an approach to calculate fluid flow in porous media with stochastic distribution of fracture.

For the case of random orientation of fractures (overall isotropy), the average value of Self-consistent scheme consists in replacing the matrix permeability k in the fracture contribution k f (second term in the right hand side of equation ( 39)) by the unknown effective permeability k = k sc . This leads the following second degree equation for k sc . Indeed, the dilute Mori-Tanaka's estimate [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF] corresponds to the first order development of Equation ( 40) with respect to ρ.

n n ⊗ is 1 3 δ δ δ δ and
In order to unify the results for other fracture shapes, it is useful to use dimensionless variables. We recall the definition of the excluded volume V ex to define the dimensionless density. The excluded volume for a fracture is defined as the volume into which the center of another fractures may not enter if there is no fracture intersection. Therefore, the dimensionless density ρ that is defined as the number of fractures per excluded volume reads [START_REF] Adler | Fractures and Fracture Networks[END_REF]:

ex V ρ = ρ (41)
where, V ex =π 2 R 3 for disc fracture of radius R [START_REF] Charlaix | A criterion for percolation threshold in a random array of plates[END_REF]. In order to study the variation of effective permeability versus ρ , another dimensionless variable is introduced as follows: which is exactly the singular value of density obtained by Pouya and Vu [START_REF] Pouya | Flow around a disc-shaped crack inclusion in infinite matrix. Application to permeability upscaling of fractured porous media[END_REF] for the case of superconductive circular fracture (i.e., c→∞). This result can have another interpretation: The limit case of λ→∞ corresponds to an impermeable matrix (k→0). For this case, the effective permeability can be different from zero beyond the critical value c ρ . According to Dormieux and Kondo [START_REF] Dormieux | Approche micromécanique du couplage perméabilité endommagement[END_REF], this critical density could be considered as the percolation threshold for the fracture network.

As it can be observed in Fig. 18, the variation of λ versus ρ , for the case k = 0 and c ρ ≥ ρ , can be approximated by a line with the following equation:

c ρ ≥ ρ ; 2 1 9 2.55 32   π = ρ -   λ   (43) Furthermore, c ρ ≥ ρ ; 2 9 2.55 32 sc k k   π = λ ρ -     (44)
For a finite value of λ , the slope of curves 1/ λ increases with ρ up to an asymptotic value. It also exits a critical value for each curve beyond which the effective permeability increases up to one order of magnitude. This critical point can be determined as the intersection of the initial and final tangents of the curve. As mentioned earlier, the effective permeability of a fractured porous medium can be evaluated by considering only a single fracture. Hence, the computation can be carried out by employing only analytical and semi-analytical elementary integrals in numerical procedure since all collocation points are located on the fracture plane. This leads to a fast and efficient numerical programme for computing equivalent permeability of a fractured material.

CONCLUSION

Numerical computation based on singular integral equations method is described to resolve the general potential solution presented by Pouya [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral Equation[END_REF] for the fluid flow in an infinite fractured porous media under a far-field pressure gradient. In this framework, a consistent mesh is only required on the surface of fractures of the studied network. Numerical solution is successfully validated by comparing numerical results with analytical solutions of two benchmark problems, namely: a single superconductive circular disc-shaped fracture in 3D and a single superconductive fracture line in 2D. A semi-analytical solution for infiltration through the fracture versus fracture conductivity, matrix permeability and fracture geometry parameters based on numerical solutions is obtained for a single fracture and different values of fracture conductivity. This is a key result to upscale the effective permeability of an infinite porous medium containing a distribution of fractures with finite conductivity.

The numerical procedure presented in this study is valid for intersecting curved fractures surface in various scales. Thus, it can be applied to model the flow in a set of curved faults under a far-field pressure condition. In particular, the developed model can be used for a more accurate simulation of hydromechanical behaviour of the "fault zones" where a dense fracture networks surround the central part of the fault [START_REF] Seyedi | A two-scale model for simulating the hydromechanical behavior of faults during CO2 geological storage operation[END_REF]. This possibility opens a way to motivate further works on the industrial application fields as petroleum, geothermic, geological CO 2 storage, etc.

It is worth reminding that only fractures with non-pressure jump between their two faces are considered in the present work. However, a fracture can act as an impermeable membrane with a pressure discontinuity across its surface. A closed-form solution for the potential discontinuity across an impervious disc embedded in an infinite matrix is proposed by Martin [START_REF] Martin | On potential flow past wrinkled discs[END_REF]. The study on 2D problem, in order to upscale the effective permeability of a porous media containing impervious crack lines, is an ongoing work and will be addressed in a future paper. The transformation from a triangular element in the global coordinate system to a reference element is described as:

0 0 . x x y y - η     =     - ξ     J or ( ) 0 .
x x ξ = -J [START_REF] Huseby | Geometry and topology of fracture systems[END_REF] with, J the matrix of transformation and (x 0 ,y 0 ) the origin of transformation. The functions thus are obtained in the global coordinate system:

( ) 0 ( ) . . k k k x x x ψ = λ + α - J (k=1,2 ,3) (48)
Interpolation functions of pressure on the elements of type 1 and 2 can be deduced from equation [START_REF] Pouya | Permeability of 3D discontinuity networks. New tensors from boundary-conditioned homogenisation[END_REF]. Their formulations in reference coordinate systems are given in equations ( 49) and (50), respectively. As a consequence, the extension of these formulations in global coordinates systems is carried out by using the following equation.

( , )

.

k k k k η     α ψ η ξ = λ +       β ξ ξ     (k=1,2 ,3) (49)
( )

2 2 2 3 2 ( , ) . 3 k k k k η     α   ψ η ξ = λ +     β -ξ     (k=1,2 ,3)

Fig. 1

 1 Fig. 1 Porous medium containing intersecting fractures In absence of point sources, continuity equation of fluid in the matrix reads:

Fig. 3

 3 Fig. 3 Masse exchange between matrix and fracture [24]: (a) at a regular point and (b) on an intersection line

2 1 r

 21 perpendicular to the fracture border.

Fig. 4

 4 Fig. 4 Three different types of element

18 )Fig. 6

 186 Fig. 6 Relative position between collocation point x c and element E (ϒ is supposed be Oxy): (a) x c in element E; (b) is on plane ϒ and not in the element E; (c) x c is not on ϒ and x cp in element E; (d) x c is not on ϒ and xcp not in element E.

Fig. 7

 7 Fig. 7 Parameters of each sub-triangle T l

Fig. 8

 8 Fig. 8 An isolated superconductive disc in an infinite porous medium

Fig. 11 Fig. 12

 1112 Fig. 11 Studied 3D band fracture and equivalent 2D configuration

Fig. 13 Fig. 15

 1315 Fig. 13 Mesh of several intersecting fractures in an infinite isotropic matrix under a far-field pressure gradient A = (1,0,0)

Fig. 16

 16 Fig. 16 Pressure field in the fracture (λ = 2)

Fig. 17

 17 Fig. 17 Numerical solution for c/ Q * (X) depending on λ, fitted by a linear function

Fig. 18

 18 Fig.18displays the variation of self-consistent estimate k sc with the fracture density for different values of fracture conductivity. In the limit case λ→∞, the self-consistent model (40) reveals a critical value

Fig. 18 Fig. 19

 1819 Fig. 18 Self-consistent scheme estimation of dimensionless variable 1/ λ plotted as a function of the dimensionless density for different values of λ

ACKNOWLEDGEMENTS

The work presented in this paper is co-funded by BRGM through an "Institut Carnot" research fund and IFSTTAR. The authors also acknowledge the fruitful discussions with Dr Pierre Sochala.

APPENDIX A

The shape functions of an arbitrary element are obtained from ones on the reference element. The equilateral triangle is selected as the reference element in the parameter space {η,ξ}(Fig. 20) with the shape functions as:

-η -ξ ψ η ξ = ;


Thomas Coisnon
email: thomas.coisnon@agrocampus-ouest.fr

Walid Oueslati
email: walid.oueslati@ncl.ac.uk

Julien Salanié

Spatial targeting of agri-environmental policy and urban development

Keywords: JEL classification : R14, R21, Q18, Q24, Q25 Land development, Urban sprawl, Leapfrog, Land rent, Monocentric model, Farming, Agri-environmental policy, Spatial targeting, Agricultural pollution

123456789A698BCDCEF9A499799

1 234564715489A56B91CD14986EABF68CBAB54713C7614B184B1AFA7C3AB5 !*D35F1DEFBDB41(56E18.8F65E1819.6E8B1565BE71 :D6.38B181A5/5E61.11B;1'230203<1 27/AE8A1'231 65FFE2E65EDB19=1#1241'241>2?41>'41>'@1 1 DF0667F1#1 5AE2E6E56EF5EDB1 8F1 8AA8F41 756838B1 .A)5EB41 3E51841 A8B81 2DB6EAA841 3DA681 3DBD68BAEB.841 51AE6.6.A841+D6EEB.8F151AE08B/EADBB838B568F416E)65181F+5E5641+D66.EDB151AE6D68"1 C87-DAF#1 65B1 8/86D+38B41 .A)5B1 F+A5-641 685+2AD141 65B1 A8B41 3DBD68BAE61 3D8641 25A3EB141 51AE0 8B/EADB38B561+D6E6741F+5E5615A18EB14151AE6.6.A561+D66.EDB"1 11 1: Il existe de manière assez répandue un soutien public destiné à la préservation d'aménités naturelles, telles que des lacs, des rivières ou des zones humides, qui s'accompagne d'efforts pour les protéger de la pollution agricole. Les politiques agri-environnementales contribuent à ces efforts en incitant les agriculteurs à adopter des pratiques respectueuses de l'environnement au sein de ces écosystèmes. Une politique agri-environnementale territorialisée favorise la préservation d'aménités naturelles, et peut ainsi influencer le choix de localisation résidentielle des ménages. L'objectif de ce papier est d'étudier l'étendue de ces effets sur une structure urbaine spatialisée. Nous utilisons un modèle de ville monocentrique dans lequel nous incluons le comportement des agriculteurs vis-à-vis d'une politique agri-environnementale. Nous constatons que l'introduction d'une politique agri-environnementale territorialisée peut mener à l'apparition de développement urbain supplémentaire, qui peut se révéler contradictoire par rapport à l'objectif même de la politique. . 5845 Widespread public support exists for the provision of natural amenities, such as lakes, rivers or wetlands, and for efforts to preserve these from agricultural pollution. Agri-environmental policies contribute to these efforts by encouraging farmers to adopt environmentally friendly practices within the vicinity of these ecosystems. A spatially targeted agri-environmental policy promotes natural amenities and may thereby affect household location decisions. The purpose of this paper is to investigate the extent of these impacts on the spatial urban structure. We extend a monocentric city model to include farmers' responses to an agri-environmental policy. Our main findings are that the implementation of a spatially targeted agri-environmental policy may lead to some additional urban development, which could conflict with the aim of the policy.

.1 1 !*D35F1DEFBDB1 1AD653+.F18.8F1 4BE/8AFE715B18AF11 *D35F"6DEFBDBD51AD653+.F0D.8F"2AE1 1 (56E18.8F65E1 1AD653+.F18.8F1 4BE/8AFE715B18AF1 8BA812DA1.A5616DBD371 4BE/8AFE71D218

Introduction

Agri-environmental policies focus on protecting and preserving natural resources and ecosystem services, although they differ throughout the world in terms of design and level of importance. Despite their international variability, agri-environmental policies (AEPs) share some common features: farmers are enrolled in environmental programs on a voluntary basis in exchange for payment [START_REF] Wunder | Taking stock: A comparative analysis of payments for environmental services programs in developed and developing countries[END_REF][START_REF] Baylis | Agri-environmental policies in the EU and United States: A comparison[END_REF].

AEPs implementation has raised several efficiency concerns. First, informational asymmetries between the contractors (i.e. farmers as 'agents' and the State as 'principal') can induce inefficiencies. These arise either due to the fact that farmers' implementation costs are unknown to the principal (adverse selection: e.g. [START_REF] Latacz-Lohmann | Auctioning conservation contracts: a theoretical analysis and an application[END_REF][START_REF] Canton | Agri-environmental schemes: Adverse selection, information structure and delegation[END_REF], or because the principal cannot control farmers' 1 actions after contract implementation (moral hazard: e.g. [START_REF] Choe | Compliance Monitoring and Agri-Environmental Policy[END_REF][START_REF] Ozanne | Moral hazard, risk aversion and compliance monitoring in agri-environmental policy[END_REF]. Secondly, as a result of budget constraints, not all of the land can be included within an AEP. This calls for targeting strategies. [START_REF] Babcock | Targeting tools for the purchase of environmental amenities[END_REF] have shown the superiority of benefit-cost targeting strategies over environmental benefits or land costs targeting strategies. When land development is frequent, land-use change should also be accounted for within an AEP spatial targeting strategy, in order to secure benefits from future land-use conversions [START_REF] Newburn | Economics and Land-Use Change in Prioritizing Private Land Conservation[END_REF][START_REF] Newburn | Habitat and open space at risk of land-use conversion: targeting strategies for land conservation[END_REF].

The efficiency of AEP targeting is highly dependent upon the correlation between costs and benefits. In this paper, we argue that this is likely to be the case in urbaninfluenced areas and that this correlation is positive, meaning that parcels of land that need the most regulations are also the most expensive ones. We propose a theoretical analysis of the interplay between urban development and agri-environmental policies.

At the rural-urban fringe, farmland is more expensive because it capitalizes the option value of future land development [START_REF] Capozza | The fundamentals of land prices and urban growth[END_REF]. This has consequences in terms of farming systems. Close to cities, farmland is substituted by non-land inputs1 and farming is more intensive [START_REF] Beckmann | Von Thünen Revisited: A Neoclassical Land Use Model[END_REF][START_REF] Heimlich | Agricultural adaptation to urbanization: Farm types in northeast metropolitan areas[END_REF][START_REF] Cavailhès | Les effets de la proximité de la ville sur les systèmes de production agricoles[END_REF]. Many ecosystem services such as open spaces, aesthetic and cultural landscapes, ground and surface water quality are by-products of farming [START_REF] Hodge | Multifunctionality in Agriculture: Evaluating the Degree of Jointness, Policy Implications. Paris: Organization for Economic Cooperation and Development. Chap. To What Extent are Environmental Externalities a Joint Product of Agriculture? Overview and Policy Implications[END_REF].

These ecosystem services may be in jeopardy near cities. Further away, farming is more extensive and agriculture produces more positive, and less negative, externalities. This spatial distribution of farming externalities also influences residential development which, in turn, modifies ecosystems. [START_REF] Irwin | Land use externalities, open space preservation and urban sprawl[END_REF] have shown that land conservation policies that aim to preserve landscape and natural amenities, have an indirect side-effect in terms of the development of surrounding parcels of land. These results were also supported by [START_REF] Roe | The effects of farmland, farmland preservation, and other neighborhood amenities on housing values and residential growth[END_REF] using a conjoint analysis, and Towe (2010) and [START_REF] Geniaux | aluation des effets des zonages environnementaux sur la croissance urbaine et l'activitricole[END_REF], respectively, studying US and European cases using a propensity score matching method. Consequently, in our approach, AEPs may affect both farming (directly) and urban development (indirectly) and the net effects of AEPs on ecosystems is undetermined.

We built a model of a monocentric city where this interplay is made explicit. Existing urban economics models assume that amenities are either exogenous [START_REF] Brueckner | Why is central Paris rich and downtown Detroit poor?: An amenity-based theory[END_REF][START_REF] Wu | The influence of public open space on urban spatial structure[END_REF][START_REF] Wu | Environmental amenities, urban sprawl, and community characteristics[END_REF] or proportional to agricultural land share [START_REF] Cavailhès | The periurban city: why to live between the suburbs and the countryside[END_REF][START_REF] Bento | Welfare Effects of Anti-Sprawl Policies in the Presence of Urban Decline[END_REF]. We modify these models by allowing the level of amenities to vary spatially with farmers' behaviour. In this endogenous setting, farming is more intensive and more polluting close to the city. We then introduce an AEP which is spatially targeted to protect a given watershed in the urban-influenced area. The implementation of the AEP depends on its adoption by farmers, which in turn depends on the opportunity cost of the land. When adopted, the AEP enables water pollution from agriculture in the watershed to be regulated. However, it also increases the attractiveness of the watershed for residential development. Inevitably, and in accordance with empirical findings by [START_REF] Hascic | Land use and watershed health in the United States[END_REF] and [START_REF] Atasoy | Estimating the effects of urban residential development on water quality using microdata[END_REF], the subsequent urban development lowers the environmental efficiency of the AEP. Thus, our model is in line with land use models and empirical evidence developed by [START_REF] Newburn | Habitat and open space at risk of land-use conversion: targeting strategies for land conservation[END_REF] and [START_REF] Langpap | Protecting watershed ecosystems through targeted local land use policies[END_REF]. We provide a rigorous framework with which to analyse these policies in an urban setting. Contrary to [START_REF] Wu | Optimal land development with endogenous environmental amenities[END_REF], who provide an insightful analysis of the dynamic of a city, our approach is static, focusing on the AEP effects. We believe that our analysis is of particular importance to European countries for two reasons. First, AEPs are poorly targeted in Europe [START_REF] Uthes | Spatial targeting of agrienvironmental measures: cost-effectiveness and distributional consequences[END_REF]. While the law offers the opportunity for Member States to target benefits and costs 2 , most European agri-environmental policy is focused on Habitat Protection Zones (Natura 2000) and Nitrate Pollution Vulnerable Zones designated by local (subnational) authorities. With the exception of a few specific cases [START_REF] Kuhfuss | Une alternative originale pour l'allocation de contrats agro-environnementaux: l'appel à projets de l'Agence de l'eau Artois-Picardie[END_REF], AEPs in Europe are not cost-benefit targeted. Secondly, European countries are densely populated and urban sprawl is an important concern [START_REF] Eea | Urban sprawl in Europe : The ignored challenge[END_REF]. The remainder of the paper is as follows. In the first section, we develop a static monocentric model of an urban area with an AEP. In the second section we present a numerical application of the model. We discuss the policy implications of our findings in the final section, and conclude by evoking some ideas for future work.

A model of location decision

This section develops a model that conforms to the basic assumptions of the monocentric city model, including an exogenously determined central business district (CBD) to which households commute for employment. Households have identical incomes and preferences, and commuting costs depend on the distance between the residential location and the CBD. Land developers have identical technologies, and the market for residential development is competitive. The model assumes a thünenian organisation of suburban agriculture, so that farmers' behaviour is influenced by the city.

The landscape is represented by a cartesian coordinate plane R 2 , with the CBD located at (0, 0) and the horizontal and vertical axes representing the West-East and North-South directions, respectively. The distance between the CBD and any residential site located at (u, v) is given by x = (u 2 + v 2 ), which together with the angle of incidence θ also determines the location of agents, where (u, v) = (x cos θ, x sin θ). The landscape is characterised by an area-featured amenity source such as a lake. All residential sites are differentiated by their proximity to this lake. Farmers located within the lake's watershed produce a residual pollution that affects the quality of amenities 2 Art. 39-4 of the European Council Regulation n o 1698/2005 provided by the lake.

The household location decision

Households make a trade-off between accessibility to the CBD and land consumption.

Each household chooses a combination of residential space q h , location (u, v), and a numeraire non-housing good s to maximise their utility subject to the budget constraint w = r(u, v)q h + s + τ x; where w is the gross household income, τ is the round-trip commuting cost per kilometre and r(u, v) is the housing rent at (u, v). Households derive utility from the amenities a (u, v) provided by the lake ecosystem (i.e. scenic lake views, fishing, recreational activities, etc).

As indicated by numerous hedonic studies on property values, the level of amenity at each location in the landscape depends on its distance from the source of the amenity.

Let a L be the maximum level of amenity provided close to the lake (u L , v L ) and z be the distance between the household's location (u, v) and the lake. The level of natural amenity at any location (u, v) is given by:

a(u, v) = 1 + a L e -ηz (1)
where η is the decreasing rate of amenity as distance from the lake increases and 1 is the level of amenity at locations that are distant from the lake.

The lake is subject to pollution produced both by farms and developed parcels of land, when they are located within the lake's watershed. The total flow of pollution that reaches the lake reduces the level of amenities that the lake ecosystem provides.

Therefore, the maximum level of lake amenities can be written as

a L = a 0 -E(W ) (2)
where a 0 is the original level of lake amenities (without pollution). The specification of the total flow of pollution E(W ) depends on the farmers' behaviour, discussed in section 2.4, and on the amount of urbanised land within the watershed.

The household utility function is assumed to be a Cobb-Douglas: U (q h , s, a(u, v)) = q β h s 1-β a(u, v) γ where 0 < β < 1 and γ > 0. The first-order conditions for the utility maximisation problem define the optimal choice of housing space and non-housing space at any location:

s * (u, v) = (1 -β)(w -τ x) (3) q * h (u, v) = β(w -τ x) r * (u, v) (4)
(3) shows that expenditure on the numeraire good decreases by (1 -β) τ with each unit of distance x from the CBD. Using (1), we can derive the houshold bid-rent function:

r * (u, v) = β β (1 -β) 1-β (w -τ x) V 1 β 1 + a L e -ηz γ β (5)
The bid-rent function (5) corresponds to a household's maximum willingness to pay for housing at location (u, v). At equilibrium, households are indifferent to where they locate because their equilibrium level of utility is the same at each location and exogenous from the perspective of an open city. In a case where there are no amenities, the bid-rent function would strictly decrease with distance to the CBD. However, with a L e -ηz > 0, households may be willing to pay more for a desirable location than for proximity to their place of work.

Developers' decision

Developers produce housing goods under constant return to scale in a competitive market, by converting agricultural land into residential use. They choose the location (u, v) and housing density d(u, v), maximising their profit function

π d (u, v) = r * (u, v)d(u, v)-r d (u, v)-c(d(u, v)), where r * (u, v) is the rent paid by households, r d (u, v)
is the land price paid by developers at the location (u, v), and c(d(u, v)) is a technical and labour cost function depending on the chosen housing density. The developer's cost function is specified as follows: c(d) = d(u, v) δ where δ > 0. Note that developers have no preferences regarding the development location. First order conditions allow us to define the optimal housing density distribution:

d * (u, v) = r * (u, v) δ 1 δ-1 (6)
Competition drives profits to zero, so that we derive the following developers' bid-price for a parcel of land located at (u, v):

r * d (u, v) = (1 -δ) r * (u, v) δ δ δ-1 (7)

The farmers' behaviour

Farmers produce a level of agricultural output per hectare y and sell it at the CBD with transport costs t per unit of distance. The production function y = f (k, (u, v)) is increasing and concave and has constant returns to scale. k (u, v) is the level of non-land inputs per unit of land used by farmers at location (u, v) and represents the intensity of farming. The prices of crops and inputs, respectively p and p k , are taken as given and the land rent at (u, v) is r a (u, v). The agricultural production function is assumed to be a Cobb-Douglas: y = Ak α with A > 0 and 0 < α < 1. In the absence of any public policy, the farmer's program is given by:

max k π a (k, (u, v)) = max k (p -tx)f (k(u, v)) -p k k(u, v) -r a (u, v) (8)
First-order conditions allow us to determine the level of non-land inputs used by farmers at distance x from the CBD:

k * (u, v) = αA(p -tx) p k 1/1-α (9) k * (u, v
) is a monotonous decreasing function of distance to the CBD. Close to the city, farms tend to be intensive, characterised by a relatively high level of non-land inputs.

Further away from the city, farms progressively become more extensive, as k falls. This spatial variation of k is consistent with many empirical observations [START_REF] Katzman | The Von Thuenen Paradigm, the Industrial-Urban Hypothesis, and the Spatial Structure of Agriculture[END_REF][START_REF] Heimlich | Agricultural adaptation to urbanization: Farm types in northeast metropolitan areas[END_REF][START_REF] Cavailhès | Les effets de la proximité de la ville sur les systèmes de production agricoles[END_REF].

Competition drives farmers' profits to zero, so that we can derive their bid-rent function at any location (u, v):

r * a (u, v) = (p -tx) 1 1-α A αA p k α 1-α (1 -α) (10)
From ( 10), we derive that ∂r * a /∂x < 0, meaning that r * a (u, v) is a decreasing function of distance. This negative bid-rent gradient is due to the Thünenian behaviour, supported by empirical evidence [START_REF] Plantinga | The effects of potential land development on agricultural land prices[END_REF][START_REF] Cavailhès | Urban influences on periurban farmland prices[END_REF][START_REF] Livanis | Urban sprawl and farmland prices[END_REF][START_REF] Wu | The Effect of the Conservation Reserve Program on Land Values[END_REF]. We assume that from a distance x, suburban agriculture stops in favour of exporting farming, preserved from any urban influence, not prone to any land-use competition, and characterised by an exogenous agricultural rent r a . Land being allocated to the highest bidder, x is reached when: r * a (x) = r a .

Lake pollution

We assume that farms located within the lake's watershed (W ) generate pollution (i.e.

sediments and pesticide run-off ). Part of this pollution ends up in the lake, and the rest is absorbed naturally as it is transported to the lake. Let φ a be the rate of pollution generated per unit of farm intensity k and λ a be the natural absorption coefficient, depending on the bio-physical characteristics of the soil (eg. texture, porosity, slope, amount of organic matter). Therefore, the net flow of agricultural pollution affecting the lake is given by:

E a (W ) = W φ a k * (u, v)e -λaz dudv ( 11
)
where z is the distance between the farm and the lake. One can reasonably assume that farms located outside the watershed have no polluting effects on the lake. In addition, as in [START_REF] Wu | Optimal land development with endogenous environmental amenities[END_REF] and [START_REF] Chen | Dynamic modeling of environmental amenity-driven migration with ecological feedbacks[END_REF], we consider that developed land also pollutes the lake. The flow of residential pollution that reaches the lake depends on the housing density of the developed area within the watershed. It is given by:

E d (W ) = W φ d d * (u, v)e -λ d z dudv (12)
where d * (u, v) is given by ( 6) and λ d represents the decreasing rate of pollution with distance from the lake.

Thus, the total flow of pollution that reaches the lake is:

E(W ) = E a (W ) + E d (W ) (13)

Implementation of a spatially targeted agri-environmental policy (AEP)

The government introduces an agri-environmental policy (AEP) aimed at reducing pollution of the lake by agriculture. The AEP only targets farmers settled within the lake's watershed. This policy consists of a voluntary contract with farmers. Contracting farmers reduce their use of non-land inputs k down to a fixed non-land inputs target level k.

In return, participating farmers are granted a uniform subsidy σ in euros per hectare.

The policy is fully described by the couple σ, k . Participating farmers' profit function is given by: max

k (p -tx)f (k(u, v)) -p k k(u, v) -r a (u, v) + σ s.t. k (u, v) ≤ k and (u, v) ∈ W (14)
There is a critical distance from which all farmers already comply with the maximum level of non-land inputs k decided by local authorities. This critical distance is given by k * ( x) = k, from which we derive:

x = 1 t p - p k αA k 1-α (15)
A rational farmer agrees to participate if the subsidy compensates for his loss of profit. The condition for a farmer to participate can be written:

(p -tx)f (k * ) -p k k * -r a ≤ (p -tx) A k α -p k k -r a + σ (16)
There is a location x ′ , from which farmers agree to participate in the program, although they have to change their practices. x ′ solves:

A αA p k α 1-α (1 -α) p -t x ′ 1/1-α -A k α p -t x ′ + p k k -σ = 0 (17)
Empirical studies confirm that the probability to participate in agricultural land preservation programs increases with the distance between the farm and the nearest city [START_REF] Lynch | Combining spatial and survey data to explain participation in agricultural land preservation programs[END_REF][START_REF] Duke | Participation in agricultural land preservation programs: parcel quality and a complex policy environment[END_REF] Figure 1: Agricultural intensity gradient and decision to adopt the policy We suppose that there are no information asymmetries between farmers and the government, and that administrative transactions are without cost. Thereby, the government has a perfect foresight of farmers' input use at any location. Figure 1 shows how the agricultural intensity is distributed across space, and which farmers, depending on their location, are likely to adopt the measure. Farmers located further than x are not granted any subsidy, although they already comply with input restrictions. Farmers located between the city and x ′ do not adopt the measures as the subsidy does not compensate them for their loss of profit. Farmers between x and x ′ within the lake watershed agree to adopt the measure and are granted the subsidy.

The spatial distribution of farming intensity is then given by:

k (u, v) =    k, if x ∈ [ x ′ , x] and (u, v) ∈ W k * , otherwise (18)
Thus by choosing k, the government determines the adoption distance x ′ and by choosing σ, sets the depth ( x -x ′ ) of the regulated area. By using (10) and ( 18), we derive the farmers' bid function:

r * a (u, v) =      (p -tx)A k α -p k k + σ, if x ∈ [ x ′ , x] and (u, v) ∈ W (p -tx) 1 1-α A αA p k α 1-α (1 -α) , otherwise (19)

Spatial equilibrium

Having derived the behavioural function for households, developers and farmers, we now characterise spatial equilibrium for land development. Six conditions combining households', farmers' and developers' decisions determine the spatial equilibrium.

1. Housing prices are bid up until no household has any incentive to move. This condition is satisfied when housing prices are represented by ( 5) since the household's bid function is the maximum willingness to pay for housing.

2. At each location the price households are willing to pay for housing equals the price developers are willing to accept for housing (r (u, v) = r * (u, v)).

3. Land prices are bid up until profits equal zero everywhere and developers are indifferent to the location of development.

4. The number of households at any location n (u, v) equals the development density divided by floor space per household : n (u, v) = d * (u, v) /q * (u, v). Land will be developed if the developer's bid price for land exceeds the agricultural rent (19).

Thus the developed area is

D = {(u, v)|r * d (u, v) ≥ r * a (u, v)}
and the total number of households (N ) in the city is given by:

D d * (u, v) q * h (u, v) dudv = N (20)
5. The boundaries of the developed area are given by a set of points (u, v), defined by: r

* d (u, v) = r * a (u, v).
6. The level of pollution that reaches the lake at equilibrium E * (W ) is given by ( 13) and depends on the developed element of the watershed, and on the proportion that is left to agricultural activity.

The lake watershed W is defined by

W = (u, v)|(u -u L ) 2 + (v -v L ) 2 ≤ z 2 W
, where (u L , v L ) denotes the coordinates of the lake and z W is the radius of the watershed. The agricultural area within the lake watershed is defined by the set W a = W -(D∩W ); that is the watershed area minus its developed part. The boundaries of the agricultural area influenced by the city are determined by the level of the exogenous agricultural rent, when r * a (u, v) = r a . We choose only to focus on cases where the targeted ecosystem is located within the area influenced by the city. In other words, we assume that z L +z W < x, meaning that the lake and its watershed are fully incorporated within the area under the influence of the city. We do not study the case where exporting independant farming is responsible for polluting the lake.

Simulation results

Parameter values

In this section, we discuss the parameter values used in our benchmark simulations. The households' income level and the proportion of their expenditure dedicated to housing were set using data from the INSEE (French National Institute of Statistics and Economic Studies) database3 . Transport costs were derived from the French Internal Revenue Service estimating commuting costs at around 0.4e/km. Finally, the equilibrium level of utility was set so that the urban rent at the CBD is plausible. V = 10, 000 leads to r * (0, 0) = 12.5 e/m 2 /year. With regard to the developer's parameter, we set δ = 4/3 as in [START_REF] Wu | Environmental amenities, urban sprawl, and community characteristics[END_REF].

Agricultural parameters such as the elasticity of production factor k and the ratio between output price and non-land input price were set from the Farm Accountancy Data Network4 . The exogenous exporting farm rent is derived from the average agricultural land price given by the Agreste (French institute of agricultural statistics)5 at around 5, 000 e/ha in France and assuming a discount rate of 0.05. Transport costs for farmers (t) and pollution parameters (φ a , φ d , λ a and λ d ) were set to ad-hoc values.

Parameter values are shown in Table 1. Households preferences for amenities (γ) and the distance between the lake and the CBD (z L ) are not specified as we discuss their impact in the next section.

Amenity location and urban development

We now present the equilibrium spatial pattern of development when no AEP is implemented. Therefore, we set σ = 0 and k = 0 for the moment.

Figure 2 depicts four different urban configurations that can be obtained depending on two key parameters : the distance between the lake and the CBD (z L ) and households preferences for amenities (γ). We present the equilibrium urban configuration for z L = (40, 50) and γ = (0.05, 0.12). The developed area can have the following patterns.

First, the lake has no influence on the city. This occurs where there is a large enough distance between the lake and the CBD (z L = 50) and preferences for amenities are low (γ = 0.05). The developed area has the characteristics of a symmetrical monocentric city (Figure 2a). Secondly, the developed shape is slightly distorted by the presence of the lake. This is the case when the lake is located nearer to the city (z L = 40). Household preferences for amenities are not sufficiently high to drive households to move right next to the lake (γ = 0.05). However, the flow of amenities provided by a more accessible lake makes the city extend in its direction (Figure 2b). Thirdly, we can observe a fragmented land development pattern. In Figure 2c, the lake is at a substantial distance from the CBD (z L = 50) but, as household preferences are higher than in case (a), some households decide to move right next to the lake, leading to the existence of two distinct developed areas. Finally, the last case is the one where we observe one unique area embracing both the CBD and the lake (Figure 2d). In this case, the lake is close enough to the city (z L = 40) and preferences are higher (γ = 0.12).

Model sensitivity to residential pollution parameter (φ d )

After highlighting the importance of the amenity's characteristics -i.e. the location of the amenity source and how it is valued by households -in the pattern of land development, we analyse the sensitivity of the model with regard to the residential pollution effect (φ d ) in the case of benchmark (c) (z L = 50 and γ = 0.12). Figure 3a depicts the size of the total developed area and of the disconnected developed area only, at different levels of impact from residential pollution. As φ d increases, land development decreases. This is due to the negative impact pollution has on the level of amenities provided by the lake. The equilibrium is reached when the development of another parcel of land would lead to a lower level of amenity, so that no more households would Note however that the level of agricultural pollution is also influenced by the residential pollution effect. A higher φ d means a smaller developed area, and therefore more agricultural land which is also likely to pollute the lake.

Impacts of a spatially targeted agri-environmental policy

Simulations of our spatial city model enable us to investigate the effects of a spatially targeted AEP on urban development. A whole set of simulations were undertaken for different combinations of (σ, k), characterising different spatially targeted AEPs implemented on our benchmark configurations depicted in Section 3.2.

In the case of benchmark (a), the lake is located at a substantial distance from the city and household preferences for amenities are low (z L = 50 and γ = 0.05), so that the lake has no influence on households or therefore on the shape of the developed area.

We observe no significant side effects of the spatially targeted AEP in terms of urban development.

In the case of benchmark (d), as the surroundings of the lake are already developed In Figure 4, we observe that, for some given (σ, k), the total developed area increases up to 1% for benchmark (b) and 0.4% for benchmark (c). The economic mechanism explaining this effect is that, for a given (σ, k), the flow of agricultural pollution reaching the lake can be significantly reduced, thereby providing a higher level of amenities.

Households may wish to pay more to be located near to the lake than they would have without the policy. An increase in household bid-rent makes developers raise their bidprice for land to the detriment of farmers. This mechanism logically leads to an increase in the developed area between the city and the lake. Obviously, whether or not the implementation of the policy has a positive impact on the size of the developed area is highly dependent upon the extent of the subsidy and on the location of the adopting area determined by local authorities via the decision of (σ, k). Therefore, the choice of policy is crucial in terms of the extent of the side effects that we identified. Furthermore, the budget aspect needs to be taken into consideration as the cost of the policy may vary depending on the choice of (σ, k).

Policy implications 4.1 Policy design and budget aspects

As indicated in the previous section, the effects of a spatially targeted AEP in terms of undesirable land development largely depend on the extent to which the policy is adopted by farmers, and therefore on the combinations of (σ, k) decided by policy authorities.

We have not taken into consideration the impacts of the policy in terms of budget until now.

Figure 5: Policy cost for different combinations of (σ, k) from benchmark (c)

Figure 5 depicts the cost of different combinations of (σ, k) implemented on benchmark (c) 6 . We observe that the policy cost varies widely. The minimum cost (0 e) is attained when none of the farmers located within the watershed have any interest in adopting the contract. This occurs when the maximum level of non-land inputs to be used ( k) is too restrictive and the subsidy (σ) is not high enough to compensate them for the reduction in profit. Conversely, it also occurs when the maximum level of non-land inputs is not restrictive enough, as this prevents farmers from adopting the AEP as they would not actually change their existing practices.

The cost increases rapidly as farmers start to adopt the contract. This corresponds to the fact that the non-land inputs limit is restrictive enough to ensure that farmers do change their practices, and that the level of subsidy compensates them for lost profits.

Obviously for a given (σ, k), the budget may also change depending on the size of the watershed and on the distance between the lake and the CBD, as farmers have varying land opportunity costs depending on their location. A larger watershed means a higher policy cost, as a greater number of farmers would have the opportunity to adopt the contract. Similarly, if the lake is at a distance where surrounding farmers are likely to be interested in adopting the contract, the anticipated budget will be higher. However, if the lake is either too close to the city, or too far from it, the budget will decrease as farmers will be unlikely to adopt the AEP as their practices will be either too intensive or the reverse.

In the case of benchmark (c), to ensure that absolutely all farmers located within the watershed agree to adopt the AEP, local authorities must implement a policy characterised by k = 1.8 and σ = 4, 500 e/ha, representing a total cost of more than 50 million euros. In reality, the level of a subsidy hardly ever exceeds 800 e/ha7 . Making all farmers adopt the AEP would be relatively expensive. Local authorities therefore need to choose who to target within the watershed itself.

In the following section, we compare the impacts of different combinations of the spatially targeted AEP, at constant budget.

Double spatial targeting of the policy

We assume that local authorities have to design their spatially targeted AEP with a fixed exogenous budget of 2.2 million euros 8 . Their perfect foresight in terms of farmers' behaviour at any location allows them to determine different combinations of (σ, k) to meet their budget. Formally, by choosing a given (σ, k), they also decide which area is likely to adopt the measure. We call this double spatial targeting. It means that authorities choose a first target, which is the watershed. In addition, their choice of (σ, k) to meet their budget will define a second target area within the watershed. In this section, we show how the impacts of the AEP may differ depending on this second spatial targeting.

We distinguish three specific cases depicted in Figure 6 In Table 2, we present the effects of the policy in terms of the pollution level and urban development side effects, from benchmark (c) and under three combinations of (σ, k) each one representing cases A, B and C respectively (see Figure 6), and all of them satisfying the fixed budget of 2.2 million euros.

With an identical budget, and depending on the second spatial targeting of authorities, we observe different effects of the policy (Table 2). We first consider the situation where residential development has no polluting effect. In case A, farmers who are likely to adopt the policy are located in the North-Western part of the watershed. In terms of pollution reduction, this policy has a positive impact as the equilibrium flow of pollution that reaches the lake decreases by nearly 14%, compared to a case where no policy is implemented. However, we observe the occurrence of an undesirable side effect in terms of development. The disconnected developed area increases by 3.41%, which, as explained before, is due to the improvement in the level of amenity induced by the policy. Conversely, in case C, the adopting area is located nearer to the city. The implementation of the policy allows farmers to resist development, as they are granted a subsidy, thereby limiting the increase in the size of the leapfrog area. In terms of pollution reduction, the policy has a lower effect (-2.60%). Indeed, as farmers are more intensive in this area than in case A, the policy has to be less restrictive if it targets this specific area.

Between than in case A (+ 2.17%).

In a situation where residential development has a polluting effect (φ d > 0), the results are similar (Table 2). However, the positive impact on the size of the total developed area and the disconnected developed area is counterbalanced by the residential polluting effect (as discussed in Section 3.3). The flow of pollution reaching the lake is higher when we introduce residential pollution (here φ d = 0.2), and the effects of the AEP in terms of pollution reduction are lower. This is explained by the additional developed area, creating more residential pollution than in the benchmark case (+ 4.50 % in case A), and therefore counterbalancing the decrease in agricultural pollution.

Discussion and conclusion

What are the lessons to be learned with regard to the design of spatially targeted agrienvironmental policies?

Spatially targeted AEPs are implemented to respond to specific environmental issues threatening local ecosystems. In cases where the ecosystem provides amenities valued by households when making decisions about their residential location, we highlight the potential side effects of a spatially targeted AEP in terms of suburban development.

As shown in the model, the lake itself has an influence on the size and shape of the developed area, depending on geographical characteristics -particularly its distance to the city centre -and on households' preferences for natural amenities.

Although designed to adress ecological issues, the introduction of a spatially targeted AEP may produce unexpected side effects in terms of urban development. Depending on the restrictiveness and level of subsidy characterising the policy, we may end up in a situation where land development is significantly accentuated, or conversely, to a situation where this effect is not observed. The occurrence of additional land development can be explained by the following economic mechanism. By reducing the pollution in the ecosystem, the implementation of the policy also favours the improvement of its quality as an amenity. More households may then be encouraged to locate further from the city and nearer to the lake to enjoy this increased level of amenity. However, this mechanism is counterbalanced by the fact that, by receiving the associated subsidy, participating farmers may have the opportunity to resist developers' bids for land. The relative weight of both mechanisms depends on which farmers adopt the policy within the target area W . More generally, authors have highlighted that basing an AEP solely on cost or benefits criteria is less efficient than a policy based on the benefit to cost ratio [START_REF] Babcock | Targeting tools for the purchase of environmental amenities[END_REF][START_REF] Uthes | Spatial targeting of agrienvironmental measures: cost-effectiveness and distributional consequences[END_REF]. Following the results of our model, we therefore propose the idea of double spatial targeting, which should consider both benefits -a decrease in the level of agricultural pollution and an improvement in the level of amenity -and costs -additional urban development and residential pollutionresulting from the implementation of an AEP.

The first target is the decision to direct the policy solely at potential polluters, i.e.

the designation of the watershed as the target area. The second spatial target involves taking into consideration the heterogeneity of farmers' behaviour. Depending on their optimal intensity level, farmers are likely either to adopt the policy or not. Thus, the choice of subsidy and restrictiveness of the policy is crucial. Within the watershed, participating farmers can be located either in an area which is likely to be developed first or in an area relatively safe from development. If the adopting area is not likely to be developed first, the AEP should achieve its main aim to decrease the flow of pollution, leading to a high level of benefit. However, improvements in the quality of amenities might lead to competition for land in favour of developers, in the area between the lake and the city, meaning a cost in terms of undesirable residential development. If the adopting area is likely to be developed first, the effects in terms of pollution reduction might not be significant, as farmers may already be more intensive in this area, and the policy should be less restrictive to convince them to adopt. This situation leads to a low benefit in terms of agricultural pollution. On the other hand, farmers may have the opportunity to increase their bid-rent, thereby counterbalancing the undesirable side effects in terms of suburban development.

Obviously, these two cases are extreme ones. There are intermediate situations where the AEP could achieve its aim without creating undesirable side effects. To meet these intermediate cases, authorities must therefore carefully design their policy, anticipating the areas where farmers would be likely to adopt, and evaluating the risks in terms of suburban development.

We also highlight that, taking into account residential pollution, the additional land development following the implementation of the policy also leads to an increase in the level of residential pollution reaching the lake. [START_REF] Hascic | Land use and watershed health in the United States[END_REF] have shown that water quality can be affected both by agricultural and residential pollution, in similar proportions. Authorities should therefore have an accurate idea of how residential land use might pollute the ecosystem they wish to protect, in order to avoid generating contradictory effects.

Further research would involve using a dynamic model in which the time process would be taken into account, reflecting the irreversibility of land conversion from agricultural to residential use and the complexity of pollution transport to the lake. Indeed, flows of agricultural or residential pollution involve complicated biological and physical processes that are not captured in our simplified model, which merely deals with interactions between the lake quality and the surrounding agriculture and residential development.

Figure 2 :

 2 Figure 2: Urban development configurations

Figure 3 :

 3 Figure 3: Sensitivity of the size of the developed area (a) and the level of pollution (b) regarding φ d . Benchmark (c).

Figure 4 :

 4 Figure 4: Variation in the size of the developed area following the introduction of a spatially targeted AEP

  . Case A denotes the situation where the choice of (σ, k) is of interest to farmers located in the North-Western part of the watershed. Case C consists of a choice of (σ, k) likely to interest farmers located in the South-Eastern part of the watershed, nearer to the city and thus more exposed to competition for residential land use. Finally, case B is an intermediate situation where farmers who are likely to adopt the contract are located in the middle of the watershed.

Figure 6 :

 6 Figure 6: Double spatial targeting

Table 1 :

 1 Parameter value and signification

	Description variable

Table 2 :

 2 these two extreme cases, we identify an intermediate situation (case B) where the flow of pollution decreases by nearly 9% and leapfrog development increases less Policy effects in terms of development and pollution at constant budget. Benchmark(c).

	Without residential pollution (φ d = 0)
	Case	Benchmark	A	B	C
	σ (e/ha)	0	300	380	720
	k	0	2.2	3.7	5.2
	Total developed area (km 2 )	2, 111.75	+ 0.26 % + 0.17 % + 0.07 %
	Disconnected developed area (km 2 )	161.5	+ 3.41 % + 2.17 % + 0.93 %
	Total flow of pollution	40.32	-13.89 % -8.88 %	-2.60 %
	Farming pollution	40.32	-13.89 % -8.88 %	-2.60 %
	Residential pollution	-	-	-	-
	With residential pollution (φ d = 0.2)
	Case	Benchmark	A	B	C
	σ (e/ha)	0	300	380	720
	k	0	2.2	3.7	5.2
	Total developed area (km 2 )	2, 072.25	+ 0.23 %	0.12 %	+ 0.05 %
	Disconnected developed area (km 2 )	126.25	+ 3.37 % + 1.58 % + 0.79 %
	Total flow of pollution	85.43	-5.00 %	-3.61 %	-1.21 %
	Farming pollution	46.28	-13.03 % -8.75 %	-3.00 %
	Residential pollution	39.15	+ 4.50 % + 2.48 % + 0.95%

Gears, chemicals, pesticides, etc.

http://www.insee.fr/fr/themes/ (accessed: May 22nd 2012)

http://ec.europa.eu/agriculture/rica/database/database_en.cfm (accessed: January 18th 2012)

http://agreste.agriculture.gouv.fr/ (accessed: August 3rd 2012)

In our simulations, we assume that 30 % of agricultural land is dedicated to non-agricultural uses (i.e. roads, farm buildings...)

In France, the detailed levels of subsidies for AEPs implemented in any region are available at http://agriculture.gouv.fr/pac-developpement-rural-feader (accessed: December 11th 2012), through the official Regional Document on Regional Development(Document Régional de Développement Rural 2007- 2013 )

 8 In the French Region Pays de la Loire, nearly 40 million euros are dedicated to the implementation of spatially targeted AEPs for the period 2010-2015 (Regional Document on Regional Development, Pays de la Loire)

"D.8F65ED51AD653+.F0D.8F"2A1