Axel Kr Öner
email: axel.kroener@math.hu-berlin.de

Eva Kr Öner
email: kroener@uni-landau.de

AND Heiko Kr Öner
email: heiko.kroener@math.uni-freiburg.de

Finite element approximation of level set motion by powers of the mean curvature

Keywords: mean curvature flow, level set equation, regularization, viscosity solution, finite elements AMS subject classifications. 35J60, 35J70, 35J75, 65L60, 35D40

In this paper we study the level set formulations of certain geometric evolution equations from a numerical point of view. Specifically, we consider the flow by powers greater than one of the mean curvature and the inverse mean curvature flow. Since the corresponding equations in level set form are quasilinear, degenerate and especially possibly singular a regularization method is used in the literature to approximate these equations to overcome the singularities of the equations. Motivated by the paper [29] which studies the finite element approximation of inverse mean curvature flow we prove error estimates for the finite element approximation of the regularized equations for the flow by powers of the mean curvature. We validate the rates with numerical examples. Additionally, the regularization error in the rotational symmetric case for both flows is analyzed numerically. All calculations are performed in the 2D case.

Introduction. Huisken and Ilmanen

 [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF]used the inverse mean curvature flow to prove the Riemannian Penrose inequality in general relativity. Later its level set formulation was extended to the flow by powers k > 1 of the mean curvature by Schulze [START_REF] Schulze | Nonlinear evolution by mean curvature and isoperimetrical inequalities[END_REF] who also proved a certain inequality using this flow. The paper [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF] arouse the interest for a numerical analysis of this special level set approach to inverse mean curvature flow which lead to the paper [START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF] by Feng, Neil and Prohl who introduced a finite element discretization for the level set formulation of inverse mean curvature flow as it appears in [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF]. They prove error estimates in the H 1 -norm and the L 2norm and confirm their rates by numerical examples. Furthermore, they focus on the aspect that their finite element method approximates the regularized equation (instead of the equation for level set inverse mean curvature flow) and present some numerical examples in which they study the corresponding regularization error.

The contribution of this paper is to study the finite element approximation of the regularized equation for level set inverse mean curvature flows as described in [START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF] for the different setting characterized by the fact that here we have flows by powers k ≥ 1 of the mean curvature, as considered in [START_REF] Schulze | Nonlinear evolution by mean curvature and isoperimetrical inequalities[END_REF]. We prove rates for the H 1 -and L 2 -error and confirm them by numerical examples. In the second part of the paper we study the regularization error in the rotational symmetric case for the flow by power k ≥ 1 of the mean curvature numerically, similarly as in [START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF]. We obtain rates within the range of the corresponding theoretical estimate from [START_REF] Kröner | Approximation rates for regularized level set power mean curvature flow[END_REF]. Moreover, similar to this estimate we observe that this rate improves when k ≥ 1 decreases. The third part of the paper deals with the regularization error for level set inverse mean curvature flow in a simplified rotational symmetric setting. Contrary to [START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF] we respect for an estimate of the regularization error the artificial boundary values. We confirm the obtained estimate in the simplified setting that the issue of artificial boundary values is omitted by a numerical example. The last section is devoted to some numerical examples for the non-symmetric case.

We give a short overview over some related publications. For the behaviour of the classical flows we refer to [START_REF] Andrews | Evolving convex curves[END_REF][START_REF] Andrews | Classification of limiting shapes for isotropic curve flows[END_REF][START_REF] Ecker | Mean Curvature Evolution of Entire Graphs[END_REF][START_REF] Gerhardt | Flow of nonconvex hypersurfaces into spheres[END_REF][START_REF] Huisken | Flow by mean curvature of convex surfaces into spheres[END_REF][START_REF] Schulze | Evolution of convex hypersurfaces by powers of the mean curvature[END_REF][START_REF] Schulze | Convexity estimates for flows by powers of the mean curvature[END_REF]. For level set formulations for mean curvature flow, see, e.g., [START_REF] Evans | Motion of level sets by mean curvature[END_REF][START_REF] Malladi | Level set methods for curvature flow, image enhancement, and shape recovery in medical images[END_REF][START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF][START_REF] Sethian | Theory, Algorithms, and Applications of Level Set Methods for Propagating Interfaces[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] and its interpretation as the value function of a deterministic two-person game see [START_REF] Kohn | A deterministic-control-based approach to motion by curvature[END_REF]. For applications in image processing of geometric PDEs we refer to [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Alvarez | Formalization and computational aspects of image analysis[END_REF][START_REF] Caselles | A geometric model for active contours in image processing[END_REF][START_REF] Chan | An Active Contour Model without Edges[END_REF]. For6 geometric flows describing the evolution of convex and non-convex curves see [START_REF] Angement | On affine heat equation for non-convex curves[END_REF][START_REF] Gage | Curve shortening makes convex curves circular[END_REF][START_REF] Gage | The heat equation shrinking convex plane curves[END_REF][START_REF] Grayson | The heat equation shrinks embedded plane curves to points[END_REF][START_REF] Mikula | Evolution of plane curves driven by a nonlinear function of curvature and anisotropy[END_REF][START_REF] Sapiro | On affine plane curve evolution[END_REF]. The approximation of geometric evolution equations with finite elements is considered in [START_REF] Barrett | On the Variational Approximation of Combined Second and Fourth Order Geometric Evolution Equations[END_REF][START_REF] Barrett | On the parametric finite element approximation of evolving hypersurfaces in R 3[END_REF][START_REF] Barrett | Parametric approximation of isotropic and anisotropic elastic flow for closed and open curves[END_REF][START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF][START_REF] Nochetto | Convergence Past Singularities for a Fully Discrete Approximation of Curvature-Driven Interfaces[END_REF][START_REF] Walkington | Algorithms for computing motion by mean curvature[END_REF][START_REF] Toulopoulos | Numerical methods for power-law diffusion problems[END_REF], by finite difference schemes in [START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature flow[END_REF][START_REF] Osher | Fronts Propagating with Curvature Dependent speed: Algorithms Based on Hamilton-Jacobi Formulation[END_REF], and by semi-Lagrangian schemes in [START_REF] Carlini | Convergence of a large time-step scheme for mean curvature motion[END_REF].

The paper is organized as follows. Section 2 introduces the setting of the level set flow by powers of the mean curvature (level set PMCF). Section 3 deals with the finite element approximation of regularized level set PMCF, proves error estimates and presents numerical examples. Section 4 presents numerically obtained rates for the regularization error of level set PMCF. Section 5 introduces the regularized level set inverse mean curvature flow formulation (level set IMCF), from [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF]. Section 6 shows theoretically and numerically obtained rates for the regularization error of regularized level set IMCF. The final Section 7 contains some numerical examples in which we simulate level set PMCF in the non-rotational symmetric case and we give a short description of the implementation used for the numerical computations presented in this paper. Finally, we give some remarks on an alternative level set formulation sometimes used in the literature for the mean curvature flow case.

Level set PMCF.

Let Ω ⊂ R n+1 be open, connected and bounded having smooth boundary ∂Ω with positive mean curvature which we consider as initial hypersurface. We call the level sets Γ t := ∂{x ∈ Ω : u(x) > t}, t ≥ 0, of the continuous function 0 ≤ u ∈ C 0 ( Ω) a level set PMCF, if u is a viscosity solution of (2.1) Section 2] for a definition viscosity solution in this case. If u is smooth in a neighborhood of x ∈ Ω with non vanishing gradient and satisfies in this neighborhood (2.1), then the level set {u = u(x)}, x ∈ Ω, is locally at x a smooth hypersurface and moves at x in the direction of its outer normal with speed H k where H is its mean curvature in x. Using elliptic regularization of level set PMCF we obtain the equation

     div Du |Du| = - 1 |Du| 1 k in Ω, u = 0 on ∂Ω, cf. [40,
(2.2)      div Du ε ε 2 + |Du ε | 2 = -(ε 2 + |Du ε | 2 ) -1 2k
in Ω,

u ε = 0 on ∂Ω,
which has a unique smooth solution u ε for sufficiently small ε > 0, cf. [51, Section 4]; moreover, there is c 0 > 0 such that

(2.3) u ε C 1 ( Ω) ≤ c 0
uniformly in ε and (for a subsequence)

(2.4) u ε → u ∈ C 0,1 ( Ω)
in C 0 ( Ω). We call u a weak solution of (2.1), which is unique for n ≤ 6. All the above facts are proved in [START_REF] Schulze | Nonlinear evolution by mean curvature and isoperimetrical inequalities[END_REF]Section 4] under the assumption that k ≥ 1. A weak solution of (2.1) satisfies (2.1) in the viscosity sense, cf. Section [40, Section 2].

3. Discretization and error estimate for regularized level set PMCF. In this section we present a finite element discretization of the regularized equation (2.2) and prove error estimates. We will restrict us to the case that the space dimension n + 1 is 2 or 3 and that Ω is convex. The latter is only a restriction if n + 1 = 3 since ∂Ω has positive mean curvature by assumption, cf. Section 1.

W -m,p (Ω) × W m,p 0 (Ω) (F, ϕ) → F, ϕ = F ϕ ∈ R.
Let (T h , Ω h ) be a quasi-uniform triangulation of Ω with mesh size 0 < h < h 0 , h 0 sufficiently small, and V h ⊂ H 1 (Ω h ) the finite element space given by (3.2)

V h = v ∈ C 0 ( Ωh ) : v |∂Ω h = 0, v |T linear for all T ∈ T h .
In view of the convexity of Ω there holds Ω h ⊂ Ω. A function u h ∈ V h will be also considered as a function on Ω by extending it by zero in Ω \ Ω h . Then

v h ∈ H 1 (Ω).
The variational formulation of (2.2) is given by

(3.3) Ω h Du ε h , Dv h ε 2 + |Du ε h | 2 dx = Ω h (ε 2 + |Du ε h | 2 ) -1 2k v h dx for all v h ∈ V h
where we fix ε > 0 from now on and denote the finite element solution by u ε h . For formal reason we will consider boundary tetrahedrons (boundary triangles in case d = 2) to be extended to a boundary tetrahedron with one 'curved face'. Therefore we will replace a boundary element T ∈ T h (i.e. n + 1 vertexes of T lie on ∂Ω) by T = T ∪ B with (3.4)

B = {tp + (1 -t)P p | 0 ≤ t ≤ 1, p ∈ bf },
where bf is the boundary face of T , i.e. n + 1 vertexes of bf lie on ∂Ω, and P p is the unique minimizer of dist(p, •) |∂Ω . We denote the resulting triangulation by Th . This leaves the space of finite element functions we use (namely V h ) unchanged. Note, that the boundary strip Ω \ Ω h has measure O(h 2 ).

3.2. The linearized operator. We define the linear operator L ε and its dual and state some properties. Let p > 1. We define for ε > 0 and z ∈ R n (3.5)

|z| ε := f ε (z) := |z| 2 + ε 2
and denote derivatives of f ε with respect to z i by D z i f ε . We have

(3.6) D z i f ε (z) = z i |z| ε , D z i D z j f ε (z) = δ ij |z| ε - z i z j |z| 3 ε .
We define the operator Φ ε by

(3.7) Φ ε : W 1,p 0 (Ω) → W -1,p * (Ω), Φ ε (v) = -D i D i v |Dv| ε - 1 |Dv| 1 k ε ,
where 1 p + 1 p * = 1, so that (2.2) can be written as

(3.8) Φ ε (u ε ) = 0.
We denote the derivative of Φ ε in u ε by (3.9)

L ε := DΦ ε (u ε )
and have for all ϕ ∈ W 1,p 0 (Ω) that (3.10)

L ε ϕ = -D i (D z i D z j f ε (Du ε )D j ϕ) + 1 k f ε (Du ε ) -1-1 k D z j f ε (Du ε )D j ϕ =: -D i (a ij D j ϕ) + b i D i ϕ,
where we use the convention to sum over repeated indices. The coefficients a ij and b i are in C ∞ ( Ω). Note, that the estimate (2.3) is not available for higher order derivatives of u ε uniformly in ε but since ε is fixed in the present section, this does not have an effect on the following considerations. The linear operator

(3.11) L ε : W 1,p 0 (Ω) → W -1,p * (Ω)
and its adjoint operator L * ε are topological isomorphism, cf. Corollary A.2 in Section 7. From [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]Theorem 8.5.3]

we deduce for L = L ε or L = L * ε and F ∈ W -1,p * (Ω) that there is a unique solution u h ∈ V h of (3.12) Lu h , ϕ h = F ϕ h for all ϕ h ∈ V h ,
where u ∈ H 1 (Ω) is the unique solution of Lu = F and we have the estimate

(3.13) u h W 1,p (Ω) + u -u h W 1,p (Ω) ≤ c u W 1,p (Ω) . Furthermore, if F ∈ L p (Ω) we have (3.14) u -u h W 1,p (Ω) + h u -u h L p (Ω) ≤ ch 2 F L p (Ω) .
Remark 3.1. Note, that we used the assertion of [14, Theorem 8.5.3] under slightly different assumptions, namely:

(i) We assume a right-hand side

F ∈ W -1,p * (Ω) (instead F ∈ L p (Ω)).
(ii) We consider the equation on Ω (instead of a polygonal domain) and use as discretization the triple ( Th , Ω, V h ).

Error estimate.

We have the following error estimate in the W 1,p -norm.

Theorem 3.2. For every p > n + 1 and small h > 0 there exists a constant

0 < c = c( u ε W 2,2 (Ω) , p) such that (3.3) has a solution u ε h ∈ V h satisfying (3.15) u ε -u ε h W 1,p (Ω) ≤ ch.
This solution is unique in a small W 1,p -neighborhood of u ε in V h .

Proof. We set

(3.16) Bh ρ = {v h ∈ V h : u ε -v h W 1,p (Ω) ≤ ρ}, with (3.17) ρ = h λ
for an arbitrary and now fixed n+1 p < λ < 1. We will obtain u h as the unique fixed point in Bh ρ of the operator T :

V h → V h with (3.18) L ε (w h -T w h ) = Φ ε (w h ), w h ∈ V h .
We show that (i) Bh ρ = ∅, (ii) T is a contraction and (iii) T ( Bh ρ ) ⊂ Bh ρ . (i) Let I h u ε be the interpolation of u ε , i.e. the continuous piecewise linear function on Ω h which is equal to u ε at all nodes of Ω h . We extend I h u ε by zero to a function on Ω. Since (3.19)

I h u ε -u ε W 1,p (Ω) ≤ ch we have I h u ε ∈ Bh ρ for small h. (ii) Let v h , w h ∈ Bh ρ and ξ h = v h -w h then using (3.18) we conclude (3.20) L ε (T v h -T w h ) = L ε ξ h + Φ ε (w h ) -Φ ε (v h ) = (L ε -DΦ ε (v h + Θξ h ))ξ h =: F In order to estimate F W -1,p * (Ω) which leads to an estimate of T v h -T w h W 1,p (Ω) in view of Corollary A.2 we choose ψ ∈ W 1,p * 0 (Ω) with ψ W 1,p * (Ω) ≤ 1 and estimate F, ψ .
To do so we use a mean value theorem for which we need the following auxiliary estimate

(3.21) Du ε -(Dv h + ΘDξ h ) L ∞ (Ω) ≤ Du ε -DI h u ε L ∞ (Ω) + DI h u ε -Dṽ h L ∞ (Ω) ≤ ch + cρh -n+1 p ,
where ṽh = v h + Θξ h ∈ Bh ρ and where we used an inverse estimate. The resulting estimate implies

(3.22) T v h -T w h W 1,p (Ω) ≤ c(h + ρh -n+1 p ) ξ h W 1,p (Ω) ≤ 1 4 ξ h W 1,p (Ω)
for small h.

(iii) Let w h ∈ Bh ρ . We have (3.23) T w h -u ε W 1,p (Ω) ≤ T w h -T I h u ε W 1,p (Ω) + T I h u ε -I h u ε W 1,p (Ω) + I h u ε -u ε W 1,p (Ω) ≤ ρ 2 + T I h u ε -I h u ε W 1,p (Ω) + ch.
It remains to estimate the norm on the right-hand side. We have

(3.24) T I h u ε -I h u ε W 1,p (Ω) ≤ c Φ ε (I h u ε ) W -1,p * (Ω) = c Φ ε (I h u ε ) -Φ ε (u ε ) W -1,p * (Ω)
≤ ch again by a mean value theorem estimate. In view of (3.17) we have

(3.25) T ( Bh ρ ) ⊂ Bh ρ .
Employing a duality argument as in [START_REF] Feng | Error analysis of finite element approximations of the inverse mean curvature flow arising from the general relativity[END_REF] we obtain an L p -error estimate in the following theorem.

Theorem 3.3. For p > n + 1 we have

(3.26) u ε -u ε h L p (Ω) ≤ ch 2 with c = c( u ε W 2,2 (Ω) , p) > 0. Proof.
From the definitions of u ε and u ε h we get (3.27)

Ω ∇u ε |∇u ε | ε - ∇u ε h |∇u ε h | ε • ∇ϕ h dx + Ω |∇u ε | 1 k ε -|∇u ε h | 1 k ε ϕ h dx = 0 for all ϕ h ∈ V h .
This equation can be written equivalently as (3.28)

Ω (A ε h ∇e ε h ) • ∇ϕ h dx + Ω (a ε h • ∇e ε h ) ϕ h dx = 0 for all ϕ h ∈ V h .

with (3.29)

A ε h := 1 0 D 2 f ε (∇u ε + t∇(u ε h -u ε ))dt a ε h := 1 k 1 0 f ε (∇u ε + t∇(u ε h -u ε )) 1 k -1 Df ε (∇u ε + t∇(u ε h -u ε ))dt e ε h := u ε h -u ε and for later purposes we set (3.30) Āε h := D 2 f ε (∇u ε ) āε h := 1 k f ε (∇u ε ) 1 k -1 Df ε (∇u ε ). Let ϕ ∈ W
Ω |e ε h | p dx = Ω (A ε h ∇e ε h ) • ∇ϕ h dx + Ω (a ε h • ∇e ε h )ϕ h dx.
By (3.28) we have further

(3.33) Ω |e ε h | p dx = Ω (A ε h -Āε h )∇e ε h • ∇ϕ h dx + Ω ((a ε h -āε h ) • ∇e ε h ) ϕ h dx ≤ c Ω |∇e ε h | 2 |∇ϕ h |dx + c Ω |∇e ε h | 2 ϕ h dx ≤ c ϕ h W 1,p * (Ω) e ε h 2 W 1,2p (Ω) .
In view of Corollary A.2 and (3.31) we get

(3.34) ϕ h W 1,p * (Ω) ≤ c Ω |e ε h | (p-1)p * dx 1 p * = c e ε h p p * L p (Ω) .
Using (3.32) and Theorem 3.2 to estimate e ε h W 1,2p (Ω) we conclude.

The discretization errors in the H 1 -and L 2 -norm have been calculated numerically by solving Eqn. (2.1) for the case of a unit circle as initial curve with ε = 0.1 and k = 1 on meshes of various refinements (h = 0.006, ..., 0.4). The solutions have been calculated iteratively by linearizing Eqn. (3.3) in the following way:

(3.35) Ω h V j Du ε h,j , Dv h dx = Ω h f j v h dx for all v h ∈ V h
where j ∈ N is the iteration index and V j and f j have been updated in each iteration:

(3.36)

V j := 1 j = 1, γV j-1 + (1 -γ) 1 √ ε 2 +|Du ε h,j | 2 j > 1,
f j := 1 j = 1, γf j-1 + (1 -γ)(ε 2 + |Du ε h,j | 2 ) -1 2k j > 1.
The value of γ affects the convergence of the iterations and was set to 0.1. In each iteration the linearized equation has been solved using Freefem [START_REF] Hecht | New development in FreeFem++[END_REF]. These solutions where compared to a solution obtained by solving Eqn. (2.1) in radial coordinates. Note, that due to radial symmetry Eqn. (2.1) simplifies to an ODE which was solved on a one-dimensional grid of h = 0.001. Figure 1 shows the numerical discretization errors as well as the rates proven in Theorem 3.2 and Theorem 3.3. Figure 2 shows the discretization errors u ε h -u ε 0.001 and the corresponding meshes. 4. Numerical study of the regularization error for PMCF. This section shows some numerically obtained rates for the regularization error of regularized level set PMCF. Since small values of the regularization parameter ε are difficult to handle we restrict the setting in order to have very high accuracy to the rotational symmetric case. A similar case in a rotational symmetric setting which omits the important issue of artificial boundary values in the model from [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF] was considered in [29, Figure 5] for regularized level set IMCF.

We state the following theoretical regularization error estimate.

Theorem 4.1. Let u be solution of (2.1) and u ε of (2.2). Then for λ > 2k the regularization error of the regularized level set PMCF satisfies Fig. 2. Top: mesh with decreasing meshsize h. Bottom:

|u -u ε | C 0 (Ω) = O(ε 1 λ ).
u ε h -u ε 0.001
Proof. See [START_REF] Kröner | Approximation rates for regularized level set power mean curvature flow[END_REF].

We recall a well-known interpolation lemma which enables us to deduce also an approximation error in Hölder-norms. Lemma 4.2. For 0 < β < α ≤ 1 and a function v ∈ C 0 (Ω) we have

(4.2) [v] β ≤ 2 1-β α [v] β α α v C 0 (Ω) ,
where these expressions might become infinity and

(4.3) [v] α := sup x =y |v(x) -v(y)| |x -y| α .
Since u ε is uniformly bounded in the C 1 -norm, cf. (2.3), we can use Lemma 4.2 to conclude

(4.4) u -u ε C 0,β (Ω) ≤ c(β)ε λ(1-β)
for every 0 < β < 1 and 0 < λ < 1 2k . Our rotational symmetric example for which we study the regularization error is a shrinking circle. Let the circle ∂B r0 (0) ⊂ R 2 with radius r 0 > 0 serve as our initial curve so that the solution u of (2.1) is given as

(4.5) u(r) = r k+1 0 -r k+1 k + 1 ,
where r ≥ 0 denotes the radius variable in polar coordinates in R 2 with center in 0.

The regularized equation reduces to the following one-dimensional equation when formulated in radial coordinates (4.6)

           1 r d dr   r d dr u ε ε 2 + d dr u ε 2   = -ε 2 + d dr u ε 2 -1 2k in [0, r 0 ], d dr u ε (0) = 0 and u ε (r 0 ) = 0.
We solve this one-dimensional boundary value problem by a Newton-Algorithm combined with Thomas Algorithm as direct solver. In Figure 3 , that means indeed that the rate depends on the power k -the larger k the smaller the rate. However, this is a little bit better than 1/(2k).

Figure 4 shows the corresponding solutions of the regularized equations.

5. Level set IMCF. We recall from [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF] the following facts. Let M ⊂ R n+1 be a closed, embedded hypersurface. A classical solution of the inverse mean curvature flow is a smooth family x : M × [0, T ] → R n+1 of hypersurfaces M t := x(M, t) satisfying the parabolic evolution equation (5.1)

∂x ∂t = ν H , x ∈ M t , 0 ≤ t ≤ T,
where H, assumed to be positive, is the mean curvature of M t at the point x and ν is the outward unit normal. If the flow is given by the level sets of a function u : R n+1 → R via (5.2)

E t := {x ∈ R n+1 : u(x) < t}, M t := ∂E t ,
then, wherever u is smooth with ∇u = 0, equation (5.1) is equivalent to set equation is given by (5.5)

           E ε u ε := div ∇u ε |∇u ε | 2 + ε 2 -|∇u ε | 2 + ε 2 = 0 in Ω L , u ε = 0 on ∂E 0 , u ε = L -2 on ∂F L
where L, ε > 0 and in case there exists a (hence unique) solution u ε we will denote it by u ε,L as well. From [39, Lemma 3.4] we know the following existence result.

Lemma 5.1. For every L > 0 there is ε(L) > 0 such that for 0 < ε ≤ ε(L) a smooth solution u ε of (5.5) exists.

Furthermore, [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF]Example 2.3] shows that one expects that L is at most c ε(L) , or equivalently,

(5.6) ε(L) ≤ c L .
From [39, page 365] we recall the definition of a weak solution of (5.3).

Definition 5.2. (i) Let Ω ⊂ R n+1 be an open set then u ∈ C 0,1 (Ω) is a weak solution of (5.3) on Ω if

(5.7) J K u (u) ≤ J K u (v)
where

(5.8) J K u (v) := K |∇v| + v|∇u|dx for all v ∈ C 0,1 (Ω) with {v = u} ⊂⊂ Ω and compact {v = u} ⊂ K ⊂ R n+1 .
(ii) u is a weak solution of (5.3) with initial condition E 0 if

(5.9) u satisfies (i) with Ω = R n+1 \ Ē0 and E 0 = {u < 0}.

We state the following existence theorem which holds due to [39, Theorem 3.1] and remark that the latter contains further information like a gradient estimate.

Theorem 5.3. For every open and bounded set Ω ⊂ R n+1 there is a weak solution of (5.3) with initial condition E 0 which is unique on R n+1 \ E 0 .

Furthermore, from the proof of [39, Theorem 3.1] we conclude the following. There exist R L i → ∞, 0 < ε i → 0, solutions u i = u εi,Li of (5.5) (i.e. with ε = ε i and L = L i ) and u ∈ C 0,1 (R n+1 \ E 0 ) so that (5.10) u i → u locally uniformly on R n+1 \ E 0 and u is a weak solution of (5.3). Furthermore, the proof of [START_REF] Huisken | The inverse mean curvature flow and the Riemannian Penrose inequality[END_REF]Theorem 3.1] shows that we may assume that ε i is small compared to 1 Li . 6. Barriers and regularization error for rotational symmetric level set IMCF. The regularized equation (5.5) depends on the triple (ε, Ω L , L) where Ω L and L are coupled explicitly. We derive in a first step upper and lower bounds for the solution of (5.5) depending on the data. Then we simplify the setting by assuming that E 0 is a ball and that the boundary values on ∂Ω L coincide with (exact) IMCF starting from ∂E 0 . Our general barriers imply in this special case that the regularization error is of order ε 2 which we confirm with a numerical example. We remark that in general the 'regularization error' is mainly dominated by the artificial boundary values. We think that it is interesting to study the regularization error for more general E 0 .

Let positive L and ε be given so that problem (5.5) has a solution ū = u ε, L. As stated above we want to estimate the regularization error which means here to estimate (6.1) |u -ū| C 0 ( Ωl )

for some fixed 0 < l < L, where u is a weak solution of (5.3). Since the boundary values on ∂F L in (5.5) are rather artificial, (6.1) can only be expected to be small for 0 < l L. The idea to derive an estimate for (6.1) is as follows. Let δ > 0 be arbitrary then there exists i = i(δ) ∈ N so that

(6.2) |u -u εi,Li | C 0 (Ω l ) ≤ δ and ε i ε and L i L.
We then construct barriers for (5.5) for general ε > 0 and L > 0 for which (5.5) has a solution and deduce from these bounds in the cases (ε, L) = (ε, L) and (ε, L) = (ε i , L i ) bounds for ū and u εi,Li , respectively, which imply an estimate for ū -u εi,Li . We assume that ε > 0 and L > 0 are so that (5.5) has a solution. We use the ansatz ϕ(v) where ϕ ∈ C ∞ (R) will be chosen appropriately. We have

(6.3) D i v(x) = (n -1) x i |x| 2 , D i D j v(x) = (n -1) δ ij |x| 2 -2 x i x j |x| 4 and (6.4) ∆v(x) = (n -1)(n -2) |x| 2 .
Setting r = |x| and w = (|ϕ

| 2 |Dv| 2 + ε 2 ) 1 2
we have (6.5)

E ε ϕ(v) = div ϕ Dv (|ϕ | 2 |Dv| 2 + ε 2 ) 1 2 -|ϕ | 2 |Dv| 2 + ε 2 1 2 = w -3 (ϕ |Dv| 2 w 2 + ϕ ∆vw 2 -|ϕ | 2 ϕ |Dv| 4 -(ϕ ) 3 D i vD j vD i D j v -w 4 ). Using w = ( (n-1) 2 r 2 |ϕ | 2 + ε 2 ) 1 2
we obtain further (6.6)

E ε ϕ(v) = w -3 ϕ (n -1) 2 r 2 w 2 + ϕ (n -1)(n -2) r 2 w 2 -|ϕ | 2 ϕ (n -1) 4 r 4 + (ϕ ) 3 (n -1) 3 r 4 -w 4 = w -3 (ϕ ) 3 r 4 (n -1) 4 + ε 2 r 2 (ϕ (n -1) 2 + ϕ (n -1)(n -2)) - (n -1) 4 r 4 |ϕ | 4 -2ε 2 (n -1) 2 r 2 |ϕ | 2 -ε 4 .
assumption 6.1. We assume

(6.7) ε = αL -1 ,
with 0 < α < 1 2 . In view of (5.6) this assumption is a not too restrictive and will be assumed in the following. Our aim is that the leading term of the right-hand side of (6.5) is

(6.8) (n -1) 4 r 4 (ϕ ) 3 (1 -ϕ )
and this is enforced by exploiting the fact that all remaining terms (6.9)

ε 2 r 2 (ϕ (n -1) 2 + ϕ (n -1)(n -2) -2(n -1) 2 |ϕ | 2 -ε 2 r 2 )
have the factor ε 2 . We will choose ϕ with |ϕ | < 1 and 0 < |1 -ϕ | = δ < 1 so that the term (6.8) becomes leading if (6.10)

|(ϕ ) 3 (1 -ϕ )| ≥ (1 -δ) 3 δ ≥ 12α 2 (n -1) -2 ≥ α 2 ((n -1) -2 + (1 + δ)(n -1) -2 + 2(1 + δ) 2 (n -1) -2 + 1 8 (n -1) -4 ).
We assume that ∂E 0 ⊂ B r2 (0) \ B r1 (0) with 0 < r 1 < r 2 . Our ansatz for the upper barrier is ϕ 1 (v) and for the lower barrier ϕ 2 (v) with ϕ i , i = 1, 2, a linear function with slope 1 + (-1) i+1 δ which lies above (case i = 1)/below (case i = 2) and touches the line segment which connects the points ((n -1) log r i , 0) and (L, L -2) and δ satisfies (6.10). From a comparison principle we know that these barriers are bounds for u ε from above and below in Ω L , furthermore, these bounds are obtained explicitly for given data L, α as follows. For the upper bound we obtain (6.11)

ϕ 1 (x) = s 1 (x -(n -1) log r 1 ) if s 1 > L-2 L-(n-1) log r1 , s 1 (x -L) + L -2 else,
where x ∈ R. Analogously, we have for the lower bound (6.12)

ϕ 2 (x) = s 2 (x -(n -1) log r 2 ) if s 2 < L-2 L-(n-1) log r2 , s 2 (x -L) + L -2 else.
These barriers provide good estimates in Ω l , 0 < l L, if L is large, α small and r i < 1 are both close to 1, then especially the initial hypersurface has 'small oscillation'. Note, that since the boundary values L -2 of u ε on ∂F L are rather artificial we also expect good estimates only in Ω l (and not in Ω L ). We recall that applying this construction of barriers for the pairs (ε, L) and (ε i , L i ) leads to barriers φ1 (v), φ2 (v), ϕ i 1 (v), and ϕ i 2 (v) with obvious notation and that we have then the inequality (6.13)

ϕ i 2 (v) -φ1 (v) ≤ u -u ε ≤ ϕ i 1 (v) -φ2 (v)
which is an estimate of the regularization error.

We will exploit this way to construct barriers in the special situation when E 0 is a ball with center in the origin and radius r and when we choose as boundary values on ∂Ω L those from (exact) level set IMCF. Then the requirement for the barriers is that ϕ i , i = 1, 2, is a linear function with slope 1 + (-1) i+1 δ which lies above (case i = 1)/below (case i = 2) and touches the line segment which connects the points ((n -1) log r, 0) and (L, L) where L ∈ R suitable (i.e. equal to the 'arrival time' of IMCF at the boundary of B L (0)) and δ satisfies (6.10). It is clear that this yields a regularization error which is 'purely due to ε' and which is given by ε 2 . To see this note that the constant δ involved in the definition of the slopes of the barriers is of size δ ≈ α 2 ≈ ε 2 .

In the corresponding implementation we solved the equivalent 1D-problem on the interval [0, log L] given by (6.14) (ϕ ) 3 e 4x + ε 2 e 2x ϕ -

1 e 4x |ϕ | 4 -2ε 2 1 e 2x |ϕ | 2 -ε 4 = 0 with boundary values ϕ(0) = 0, ϕ(log L) = log L which is solved by ϕ(x) = x if ε = 0.
This is illustrated in Figure 5. Note that [29, Table 5] differs by the obtained rate O(ε) instead of O(ε 2 ) as in our case. This might be due to the fact that the example in [29, Table 5] simulates IMCF backward in time, i.e. the boundary of the domain is the hypersurface which is reached by the IMCF starting from level sets in the interior of the domain. Then the solution has stationary points. We study the forward problem, i.e. we calculate what happens with a given initial hypersurface in the future and hence we don't have stationary points in the rotational symmetric case.

7. Simulations and further remarks. In this section we present some simulations in the non-rotational symmetric case for PMCF. A short description of the implementation used for the numerical examples is presented, and finally we provide an alternative level set formulation. 7.1. Simulations in a non-rotational symmetric case for PMCF. The phenomenon of becoming round can be measured by the isoperimetrical deficit (7.1)

l(t) 2 -4πa(t),
where l(t) denotes the length of the curve and a(t) the enclosed area at time t.

According to theoretical results in [START_REF] Schulze | Nonlinear evolution by mean curvature and isoperimetrical inequalities[END_REF] we confirm the monotonicity of this deficit during the evolution in the special case of the ellipse as initial curve, see Figure 6. Furthermore, we see that with increasing k and ε = 0.05 the curves transform faster into a circle except for k = 0.5, see Figure 8, where the 'point' is reached quite well).

In Figure 4 and Figure 7 we see when comparing the exact solutions for the circle for different values of k and the approximate solutions for the ellipse with ε = 0.15 for different values of k, respectively, that the flow reaches the singularity for larger k earlier.

In our previous Figure 4 we plot a section (along the long and short half axes of the initial curve) of the solution u ε in the case of the circle and in Figure 7 in the case of an ellipse as initial curve for different values of ε.

Figure 8 shows level sets of u 0.1 for the case of the ellipse as initial curve and different values of k. We remark that our theory covers only the case k ≥ 1 but there is a well-defined and well-known behavior for the flow of convex curves with speeds given by general positive powers of the curvature, see [START_REF] Andrews | Classification of limiting shapes for isotropic curve flows[END_REF]. Our observations are as follows. For k = 0.5 we see for ε = 0.1 a quite good approximation of the phenomenon of shrinking to a 'round point' and further lessening of ε does not show significant improvements. For all k the inner level line for ε = 0.1 seems to be already 'round' while for k = 2 this seems to be far from a 'point'.

On the implementation. To calculate the finite element approximation u ε

h of u ε we used a discretization with unstructured grids, see Figure 9. These were generated by the mesh generator Gmsh, see [START_REF] Geuzaine | Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. We solved the non-linear equation (3.3) with a Newton method which uses a bi-conjugate gradient stabilized solver (BiCGSTAB) and SSOR preconditioning. For the implementation we used PDELab, a discretization module for solving PDEs which depends on the Distributed and Unified Numerics Environment (DUNE). As further references concerning PDELab we refer to [48,[START_REF] Bastian | Generic implementation of finite element methods in the distributed and unified numerics environment (dune)[END_REF], information about DUNE can be found in [START_REF] Blatt | The Iterative Solver Template Library[END_REF][START_REF] Bastian | A generic grid interface for parallel and adaptive scientific computing. Part ii: Implementation and tests in DUNE[END_REF][START_REF] Bastian | A generic grid interface for parallel and adaptive scientific computing. part i: abstract framework[END_REF][START_REF] Dune | Distributed and unified numerics environment[END_REF]. In order to get solutions for small ε we used a warm-start, i.e. we decreased ε stepwise to the desired small value and performed on each stage a calculation with the solution for the previous ε as initial value.

7.3. Alternative level set formulation. For completeness we mention an alternative formulation of the motion by powers of the mean curvature which uses a level set formulation with a level set function which depends on the time. Let M 0 ⊂ R n+1 be a given initial hypersurface and u 0 : R n+1 → R a continuous function such that (7.2)

M 0 = {x ∈ R n+1 : u 0 (x) = 0}. Let u : [0, ∞) × R n+1 → R be the unique viscosity solution of (7.3) d dt u =|Du| div Du |Du| k in R n+1 × (0, ∞) with u(0, •) = u 0 in R n+1
, we call the family of the (7.4)

M (t) = {x ∈ R n+1 : u(t, x) = 0}, t > 0,
a (time dependent) level set PMCF. Equation (7.3) is a fully nonlinear, degenerate and possibly singular (if Du = 0) parabolic equation. In the case k > 1 the elliptic main part of (7.3) is not in divergence form and fully nonlinear in the second spatial derivatives which is of disadvantage having our finite element approach in mind. Furthermore, it is higher dimensional than the equation we used. Nevertheless, this formulation is quite common in the literature in the cases 0 < k ≤ 1 and in general also available when the speed is not necessarily positive. We give a short overview.

Existence and uniqueness of a solution for (7.3) is proved in [START_REF] Chen | Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations[END_REF][START_REF] Chen | Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations[END_REF][START_REF] Evans | Motion of level sets by mean curvature[END_REF] in the case k = 1. In [START_REF] Mitake | On convergence rates for solutions of approximate mean curvature equations[END_REF] equation (7.3) in case 0 < k ≤ 1 is approximated by a family of regularized equations and rates of convergence of the corresponding solutions are obtained.

Concerning regularizations of equations we also refer to [START_REF] Bellettini | Elliptic approximations pf prescribed mean curvature surfaces in Finsler geometry[END_REF]. The time dependent formulation (7.3) in the case k = 1 3 , i.e. the affine curvature equation, is used for image processing, cf. [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Giga | Surface evolution equations. A level set approach[END_REF]. In the case k = 1, i.e. mean curvature flow, equation (7.3) has been studied intensively analytically and numerically, cf., e.g., [START_REF] Carlini | Convergence of a large time-step scheme for mean curvature motion[END_REF][START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature flow[END_REF][START_REF] Deckelnick | Convergence of a finite element method for non-parametric mean curvature flow[END_REF][START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF][START_REF] Kohn | A deterministic-control-based approach to motion by curvature[END_REF]. We want to point out the paper [START_REF] Deckelnick | Error bounds for a difference scheme approximating viscosity solutions of mean curvature flow[END_REF] by Deckelnick, where the solution u ε of a regularized version of (7.3) is approximated by a finite difference scheme which was originally proposed by Crandall and Lions [START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature flow[END_REF]. In Deckelnick's paper rates for the convergence of the discrete solution to the solution u of the (not regularized) level set equation are proved. The total error consists of a regularization error of the form (7.5) u -u ε L ∞ (Ω) ≤ c α ε α with α ∈ (0, 1 2 ) arbitrary and c α a positive constant, see [START_REF] Deckelnick | Error bounds for a difference scheme approximating viscosity solutions of mean curvature flow[END_REF]Theorem 1.2] for details, and a discretization error which is a polynomial expression in the numerical parameter and the reciprocal regularization parameter. Furthermore, the value for the convergence order of the discretization error (and hence for the total approximation error) is very low; the main point here is that this rate is of polynomial order. A calculation showed that in the setting of our paper the constants involved in the L 2 -and H 1error estimate for the discretization error depend at most exponentially on 1 ε . Such an exponential dependence is not unusual as can be seen in the paper [START_REF] Deckelnick | Convergence of numerical schemes for the approximation of level set solutions to mean curvature flow[END_REF]. There the viscosity solution u of (7.3) is approximated by a solution u ε of the regularized equation and then the regularized equation is approximated by a solution u ε,h of a semi-discrete problem. The regularization error is again of the form (7.5) but the error u ε -u ε,h measured in a certain energy norm, cf. [START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF]Theorem 6.4], is only of order c ε h, where, and this is the important point, the constant c ε depends exponentially on 1 ε . Numerical tests from that reference, however, suggest that the resulting bound overestimates the error. In the special case of two dimensions, i.e. the moving hypersurfaces are curves, Deckelnick and Dziuk [START_REF] Deckelnick | Convergence of numerical schemes for the approximation of level set solutions to mean curvature flow[END_REF] prove L ∞ -convergence (without

3. 1 . 2 0(

 12 Discretization. We fix some standard notation concerning finite elements. We denote the Euclidean norm of R n by | • |. For an open subset Ω of R n and m ∈ N * , p ≥ 1, we denote the corresponding Sobolev spaces by W m,p (Ω), W m,p 0 (Ω), H m (Ω) = W m,2 (Ω) and H m 0 (Ω) = W m,Ω). The dual spaces are denoted by W -m,p (Ω) = W m,p 0 (Ω) * and the dual pairing by (3.1)

Fig. 1 .

 1 Fig. 1. Discretization error for the unit circle as initial curve with k = 1, ε = 0.1.

  (a) the error u ε -u L ∞ (Ω h )

Fig. 3 .

 3 Fig. 3. (a) regularization error in case of a circle as initial curve for k = 1, 1.2, 1.5, 2, 3; (b) rate of regularization error as function of k.

Fig. 4 .

 4 Fig. 4. Radial solution for a circle as initial curve for (a) k = 1, (b) k = 1.5, (c) k = 2 and (d) k = 3.

2 Fig. 5 .Fig. 6 .

 256 Fig. 5. Regularization error in case of a circle as initial curve.

Fig. 7 .

 7 Fig.7. Solution for ellipse as initial curve. Picture 1-3: Section in direction of the long half axis of the initial curve for k = 1, 1.5, 2; Picture 4-6: the same for the short half axis.

Fig. 8 .

 8 Fig. 8. Solution for ellipse for ε = 0.1.

Fig. 9 .

 9 Fig. 9. Mesh for the discretization with size h = 0.15 for ellipse with half axes 1 and 2.

Acknowledgment. Parts of the work were carried out while the third author benefited from a Weierstrass Stipendium of the Weierstrass Institute-Berlin (WIAS).

rates) of the discrete solution provided h = h(ε) sufficiently small, where 'sufficiently small' is not given by an explicit formula or polynomial dependence.

Appendix A. Some auxiliary observations. Since L ε : H 1 0 (Ω) → H -1 (Ω) is a topological isomorphism by classical L 2 -theory this also holds for L * ε : H 1 0 (Ω) → H -1 (Ω). We define the to L ε associated uniformly, elliptic regular Dirichlet bilinear form of order 1 by

From Fredholm's alternative, cf. [55, Theorem 10.7], we deduce that for every F ∈ W -1,p * (Ω) the equation

3) has a unique solution.

Lemma A.1. dim N p * = dim N p = 0.

Proof. Let v ∈ N p * . From [55, Theorem 7.6] we get v ∈ W 1,p 0 (Ω) for all 1 < p < ∞, especially for p = 2. Since we know from L 2 -theory that (A.3) has a unique solution u ∈ W 1,2 0 (Ω) if p = 2 and F = 0 we deduce that v = 0. Analogously we obtain the remaining claim.

By bounded inverse theorem we conclude the following result.

Corollary A.2. L ε , L * ε are topological isomorphisms.


Cécile H Albert
email: cecile.albert@imbe.fr

Bronwyn Rayfield

Maria Dumitru

Andrew Gonzalez

Applying network theory to

Keywords: biodiversity, corridors, habitat fragmentation, metapopulation, graph theory, species distribution models, Zonation

come

INTRODUCTION

Conserving connectivity [START_REF] Crooks | Connectivity conservation[END_REF] -the capacity of species to move across landscapes -is a preferred strategy for mitigating biodiversity loss due to land-use change and habitat fragmentation and for fostering the migration, natural adaptation, and increased persistence of species under climate change [START_REF] Heller | Biodiversity management in the face of climate change: a review of 22 years of recommendations[END_REF]. In regions that have already been substantially transformed by human activities, such as metropolitan areas [START_REF] Irwin | The evolution of urban sprawl: evidence of spatial heterogeneity and increasing land fragmentation[END_REF], conserving biodiversity given uncertain scenarios of climate and land-use change requires adopting a network approach that manages remnant habitat for connectivity [START_REF] Schmitz | Conserving biodiversity: practical guidance about climate change adaptation approaches in support of land-use planning[END_REF].

Incorporating connectivity objectives and climate and land-use change uncertainties into the design of protected-area networks is a major challenge [START_REF] Mumby | Reserve design for uncertain responses of coral reefs to climate change[END_REF]) that requires new methods to prioritize areas that simultaneously satisfy multiple conservation criteria [START_REF] Zetterberg | Making graph theory operational for landscape ecological assessments, planning, and design[END_REF]. Habitat networks should account for the needs of multiple species with contrasting life histories and movement ecologies, rather than the needs of a single or few umbrella or flagship species [START_REF] Nicholson | A new method for conservation planning for the persistence of multiple species[END_REF]. The network's structure should allow movement to occur at multiple spatial scales because the contributions of habitat patches to connectivity vary across species and movement types [START_REF] Minor | A multiscale network analysis of protected-area connecitivty for mammals in the United States[END_REF][START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]. Habitat networks should also account for uncertainties in climate-driven range shifts because the contribution of habitat patches to connectivity may change as species shift their ranges [START_REF] Beier | Toward best practices for developing regional connectivity maps[END_REF][START_REF] Kujala | Conservation Planning with Uncertain Climate Change Projections[END_REF].

Finally, the network's connectivity must be robust to ongoing land-use change that can alter both habitat patches and the surrounding matrix [START_REF] Coulter Riordan | Land use compounds habitat losses under projected climate change in a threatened California ecosystem[END_REF].

The integration of multiple objectives can be achieved by combining new approaches to connectivity -measurement and mapping of contributions of habitat patches to multiple connectivity criteria [START_REF] Carroll | Use of linkage mapping and centrality analysis across habitat gradients to conserve connectivity of gray wolf populations in Western North America[END_REF]) -with systematic conservation planning tools designed for multipurpose spatial prioritization [START_REF] Moilanen | Prioritizing multiple-use landscapes for conservation: methods for large multi-species planning problems[END_REF]). These prioritization tools can balance (sometimes conflicting) needs in terms of habitat quality and connectivity for different species under distinct climate and land-use projections (e.g., [START_REF] Magris | Integrating multiple species connectivity and habitat quality into conservation planning for coral reefs[END_REF][START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]. Graph (or network) theory, in which landscapes are conceptualized as networks of weighted nodes (habitat patches with different qualities) connected by weighted links (potential movement based on properties of the intervening matrix) is a powerful tool to model habitat connectivity with available data [START_REF] Calabrese | A Comparison-Shopper's Guide to Connectivity Metrics[END_REF].

Worldwide, governments are requesting research to help plan regional-scale habitat connectivity that is robust to both land-use and climate change [START_REF] Beier | Toward best practices for developing regional connectivity maps[END_REF][START_REF] Brodie | Connecting science policy and implementation for landscape scale habitat connecitivty[END_REF]. There is thus a need for a general framework to support the design of multispecies global-change-proof habitat networks. To address this challenge, we developed an integrative framework that prioritizes areas for conservation based on maintaining regional habitat quality and connectivity for a broad set of species under climate and land-use change (Fig. 1).

We focused on 2 dimensions of connectivity: short-range connectivity that alongside habitat quality facilitates the long-term persistence of multiple species within the habitat network and long-range connectivity that promotes the maintenance of seasonal and climate-driven migrations across the habitat network [START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]. This spatial prioritization can guide the implementation of protected-area networks for conservation of biodiversity and ecosystem services. Although our approach is general and could be applied to most highly fragmented regions, we used a specific study system to demonstrate its application.

We received a direct request from Quebec government to apply our framework to guide the selection of protected areas for the southwestern region of the St Lawrence lowlands around Greater Montreal, Canada (~27,500 km 2 ) (Supporting Information). This highly fragmented periurban territory is undergoing rapid sprawl of low density urban areas, but there is strong political will and commitment from diverse stakeholders to conserve habitat connectivity within and across the region to facilitate northward climate-driven range shifts expected between the Appalachian to the Laurentian Mountains [START_REF] Mitchell | The Montérégie Connection Project : Linking landscapes, biodiversity, and ecosystem services to improve decision-making[END_REF][START_REF] Dupras | The impacts of urban sprawl on ecological connectivity in the Montreal Metropolitan Region[END_REF]. This makes it a particularly suitable test situation for our framework.

We applied our framework to a broad set of species to provide a spatial conservation prioritization of the region's natural forest habitat. We compared prioritizations based on four sets of conservation criteria to show how balancing species-specific requirements for habitat quality, short-and long-range connectivity, and climate suitability modified conservation priorities. We derived different conservation scenarios from these prioritization schemes and investigated their effectiveness, relative to a business-as-usual scenario of land-use change, at maintaining habitat connectivity into the future with spatially explicit dynamic land-use simulations.

METHODS

We used 6 steps structured within a loop to identify spatial conservation priorities (Fig. 1).

First, we identified a set of vertebrate species that represented the range of responses to habitat fragmentation and climate change present in the region based on their life-history and movement traits. Second, we mapped habitat quality and parameterized a resistance surface for each species. Third, we derived habitat networks from habitat quality maps and estimated species dispersal abilities to quantify the contribution of habitat pixels to short-and longrange connectivity for each species. Fourth, we assessed changes in species-specific climate suitability under different regional climate projections with species distribution models. Fifth, we obtained spatial conservation priorities based on species-specific criteria for habitat quality, connectivity, and climate suitability with a multi-objective prioritization tool. Sixth, we established the effectiveness of different prioritization schemes into the future with spatially-explicit dynamic land-use simulations. Each step is described in detail in Supporting Information.

Selecting focal species

We selected 14 ecoprofiles in our study area; each was represented by a mammal, bird, or amphibian species [START_REF] Beier | Toward best practices for developing regional connectivity maps[END_REF]) (Table 1). We selected species that encompassed the regional diversity -regarding their connectivity and habitat needs -from a multivariate analysis based on traits known to characterize how vulnerable species are to habitat fragmentation: habitat requirements, population dynamics, and dispersal abilities (Supporting Information) [START_REF] Henle | Predictors of species sensitivity to fragmentation[END_REF]). The species also reflected the diversity of potential changes in climate suitability in temperate North America (i.e., increase, decrease or no change in area of suitable climate ) (Fig. 1 & Supporting Information), as expected with a northward range expansion for most taxa (e.g. [START_REF] Kerr | The impact of climate change on mammal diversity in Canada[END_REF].

Assessing species-specific habitat quality and resistance maps Due to the heterogeneity of available data on species occurrences and habitat in the study area, we developed habitat-quality models through a literature review (Table 1) and used raw data from multiple sources (Supporting Information). For each species, we derived a baseline habitat-quality map from a customized 8-class land-cover map at a resolution of 30 x 30 m. These baseline quality maps were then modified by a set of spatial multipliers (range 0-1) to further account for landscape composition (e.g. forest attributes) and configuration (e.g. forest edge, distance to wetlands).

Maps of habitat patches were derived from habitat-quality maps by forming groups of habitat pixels that were large enough (area>minimum patch area) and close enough (distance<gap size) to be used by a particular species (Table 1). This organism-centered definition of habitat patch reduces discrepancies between the resolution of the habitat maps (30 x 30 m) and species-specific grain perception [START_REF] Baguette | Landscape connectivity and animal behavior: functional grain as a key determinant for dispersal[END_REF] while keeping the number of habitat patches manageable (<10,000). Species-specific 5-class maps of dispersal resistance were developed in the non-habitat pixels to quantify the degree to which pixels in the matrix limit inter-patch movement relative to habitat [START_REF] Adriaensen | The application of 'least-cost'modelling as a functional landscape model[END_REF]. Resistance values were assigned based on land-cover type (e.g. intermediate in cropland, high on highways) and on the presence of linear elements (e.g. hedges). Following a scale that doubles between classes (from 1 to 32), resistance values were tailored to each species based on a literature review of their dispersal success through the different class types (Supporting Information).

Analyzing connectivity

Habitat networks were assembled by connecting habitat patches (nodes of the network) from edge to edge via least-cost paths (links of the network) through species-specific resistance surfaces [START_REF] Adriaensen | The application of 'least-cost'modelling as a functional landscape model[END_REF]. Due to the large number of nodes, we identified links between nodes with a minimum planar graph model [START_REF] Fall | A domain-specific language for models of landscape dynamics[END_REF]. Nodes and links had associated weights. Node weights included area, mean habitat quality, and mean climate suitability. Link weights included length (i.e., distance along the least-cost path), cumulative resistance, and dispersal flux. Dispersal flux between nodes i and j separated by a distance d ij (i.e., the probability that an individual in node i will disperse to node j) was calculated as a negative exponential kernel, P ij =exp(d ij x log(0.5)/D 50 ), where D 50 is the species-specific median dispersal distance. We considered 2 estimates of D 50 per species based on upper (natal dispersal distance) and lower (gap-crossing distance) limits of dispersal (Table 1).

Species-specific dispersal distances were parameterized by combining a literature review and dispersal distances estimated from species body size and life-history traits.

We used 4 connectivity metrics as input for the prioritization schemes (Table 2 &Supporting Information) [START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]. The contribution of nodes to short-range connectivity was quantified in terms of a node's importance to the total amount and quality of reachable habitat (dEC) [START_REF] Saura | Network analysis to assess landscape connectivity trends: application to European forests (1990-2000)[END_REF] and its ability to serve as a stepping stone within the network (node betweenness centrality) [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF]. The contribution of nodes to longrange connectivity was assessed based on the frequency with which they were included in the shortest path between randomly selected pairs of nodes within the Appalachian and Laurentian Mountains (modified betweenness). The contribution of each pixel to long-range connectivity was also determined based on the amount of flow through each pixel associated with moving across the landscape in multiple directions (current density) [START_REF] Pelletier | Applying Circuit Theory for Corridor Expansion and Management at Regional Scales: Tiling, Pinch Points, and Omnidirectional Connectivity[END_REF].

We used 2 network-level connectivity metrics to evaluate the effectiveness of alternative land-use planning scenarios through time (Table 2). Short-range connectivity was quantified as the total amount and quality of reachable habitat within the network, which is related to the carrying capacity of a metapopulation (EC) [START_REF] Saura | Network analysis to assess landscape connectivity trends: application to European forests (1990-2000)[END_REF]). Long-range connectivity was quantified by the conductance (easiness) associated with traversing the network and calculated as the inverse of the mean cumulative resistance of shortest paths between randomly selected pairs of nodes within the Appalachian and Laurentian Mountains (1/CR).

Modelling changes in climate suitability

We modeled changes in climate suitability for the focal species with species distribution models. After preliminary tests, we selected the five types of models with the highest predictive power and synthesized their results with a committee-averaging procedure [START_REF] Araújo | Ensemble forecasting of species distributions[END_REF]. Models were calibrated from 1971 to 2000 with continental distribution data (Supporting Information) and three uncorrelated variables (mean annual temperature, mean annual precipitations, and the ratio summer:annual precipitations) at a resolution of 10 x 10 km for all eastern North America. Models were projected for two future horizons (2025,2050) and four extreme climate-change projections selected to characterize the full range of change in precipitation and temperature for the study region in 2050 and to represent extremes of the 136 available projections derived from different climate models and emission scenarios (Supporting Information). The most likely future is thus expected to fall within them [START_REF] Beier | Toward best practices for developing regional connectivity maps[END_REF]. Climate-suitability and habitat-quality modeling were decoupled because climate-specific measures of habitat quality were not available, which precluded the investigation of interactions between these components of habitat.

Identifying spatial conservation priorities

We identified spatial conservation priorities with the conservation decision-support software Zonation version 4 [START_REF] Moilanen | Prioritizing multiple-use landscapes for conservation: methods for large multi-species planning problems[END_REF]. Zonation identifies near-optimal trade-offs among multispecies conservation criterion and provides a ranking of all pixels in the landscape from lowest to highest conservation priority. We tested four prioritization schemes with different combinations of input criteria that expanded on existing protected areas (currently 1.2% of the St. Lawrence Lowlands): habitat quality (quality based); habitat quality and connectivity (quality and connectivity based); habitat quality and climate suitability; habitat quality, connectivity, and climate suitability. Directly accounting for current and future climatesuitability layers for different climate projections in the prioritization process meant nodes that remained suitable under all climate projections were considered as relatively more important [START_REF] Kujala | Conservation Planning with Uncertain Climate Change Projections[END_REF]) (Supporting Information). For all prioritizations, we used a core-area and edge-based cell removal rule to protect core habitat for all species, increase the importance of rare features, and save computing time. All species were given equal weights.

Criteria were weighted to balance importance among habitat quality, connectivity, and climate and between short-and long-range connectivity (Table 2).

Simulating land-use change and defining land-use scenarios

To test how changes in landscape configuration affect connectivity of species habitat networks and how efficient our prioritization schemes were at maintaining short-and longrange connectivity, we simulated land-use change from 2000 to 2050. We used a model with top-down demand for new urban and agricultural areas, bottom-up constraints (e.g. soil quality), and spatial processes (e.g. diffusion) [START_REF] Verburg | Combining top-down and bottom-up dynamics in land use modeling: exploring the future of abandoned farmlands in Europe with the Dyna-CLUE model[END_REF].

We compared five contrasting land-use scenarios that differed only in the way changes were allocated spatially (Supporting Information). In the business-as-usual (BAU) scenario, The five scenarios differed only in terms of landscape configuration. All 5 were based on the same fixed rate of development: increase of 30% for urban and 4% for agriculture areas (derived from historical changes in the study region from 1990 to 2012). In our simulations, most croplands were protected from urbanization by agricultural zoning, and protected areas were fully protected from any land transformation. Habitat quality and network analyses were repeated for each species and each scenario to assess the variation in outcome through time.

Analyses were performed in R version 3.2.3 (R Development Core Team 2012); SELES version 3.4 [START_REF] Fall | A domain-specific language for models of landscape dynamics[END_REF] within the R package grainscape; ArcMap version 10.2.0.3348 (ESRI 2011);andPython version 2.7.6 (Python Software Foundation 2001-2013). The R codes are available on GitHub (https://github.com/CecileHAlbert/CB_connectivity). The BAU baseline scenario showed a loss of 12% forest cover, which translates into speciesspecific habitat loss of 10-14% by 2050. Relative to the year 2000, under BAU habitat networks were more fragmented, and there were on average 6% fewer habitat patches (range: 16% fewer -6% more) of smaller size (17% smaller -3% larger; average 5% smaller) and with more irregular edges (edge:area ratio: 3% less -15% more; average 5% more) (Supporting Information). Short-range connectivity, measured by the amount of reachable habitat (EC: 5-18% less; average 11% less), and long-range connectivity, measured by the ease of traversing the network (1/CR: 1-22% more difficult; average 12% more difficult), were reduced (Fig. 3). Habitat fragmentation in 2000 was already acute for poor dispersers (e.g. wood frog), so future land-use change caused the largest loss of connectivity for species that currently disperse well through the landscape (e.g. Pileated Woodpecker) (Fig. 3).

RESULTS

Ongoing erosion of poorly connected habitat networks

Differences among quality-based and quality and connectivity-based prioritizations

From the prioritization maps, we derived a hierarchy of conservation solutions with an increasing fraction of the area to be protected (Fig. 2). For a given amount of area protected, the quality and connectivity-based prioritization encompassed fewer and more compact blocks of land than the quality-based prioritization, including a series of large stepping-stone patches to the north that delineated the easiest paths to traverse the lowlands between the Appalachian and Laurentian Mountains (Fig. 2b). Quality-based and quality and connectivity-based prioritizations strongly differed in their highest priority areas, but the differences were weaker for lower priorities (Fig. 2a,b, Supporting Information). In both cases the protection of the top 17% of the area encompassed around 55% of the lowland forests. Little of the remaining forest could be lost without some loss of conservation value (Fig. 4, Supporting Information).

Robustness of connectivity-based prioritization to climate change

Niche modeling identified three types of potential effects on our focal species relative to climate suitability: decrease in regional climate suitability for four species, increase in suitability for one species, and relative no change in suitability for the other nine species (Fig. 1, Supporting Information).

Including present and potential future climate suitability in the quality and connectivity-based prioritization did not greatly modify the spatial conservation priorities (Fig. 2b,d, Supporting Information). Relatively few components of the top priority areas differed from when climate suitability was included (Supporting Information). Changes in top priority areas were slightly larger under quality-based prioritization and in lower priority areas (Fig. 2, Supporting Information). In both cases, a few large blocks in the northern part of the region, that offered more stable climatic conditions for coniferous-dependent species (e.g. American marten, Red-breasted Nuthatch) from 2000 to 2050 were among the top priorities when climate suitability was included in the prioritization.

Prioritizing conservation based on connectivity in a changing landscape

With the 17% quality and connectivity-based scenario, the negative effects of land-use change on connectivity were notably reduced relative to BAU (Fig. 3). Long-range connectivity, decreased slightly for all species (1/CR: 0.7 -8% less; average 3% less), whereas under BAU this decrease reached up to 22% less (Fig. 3d). Less short-range connectivity was lost in the lowlands compared with BAU (EC: 3 -15% less; average 8% less) (Fig. 3b) because the greater amount of protected-area in the lowlands caused more habitat to be lost in upland areas. This 17% quality and connectivity-based scenario led to a greater retention of the different conservation criteria: habitat quality, connectivity criteria, and even climate criteria that were not targeted in the corresponding prioritization; this was close to what would be achieved with 30% of land being protected (approximately all the forest protected) (Fig. 4).

When the top 10% of the area was protected, results differed between quality-based and quality and connectivity-based scenarios. These 2 scenarios did not significantly reduce the loss of short-range connectivity (EC: 3 -16% less and 4 -17% less respectively) and led to an intermediate loss when compared with BAU and the 17% connectivity-based scenarios (Fig. 3b). However, the quality-based scenario led to a decrease in long-range connectivity similar to BAU (1/CR: 5 -24% less), whereas the quality and connectivity-based scenario led to a significantly smaller decrease (1/CR: 0 -8% less), almost similar to the 17% quality and connectivity-based scenario (Fig. 3d). Although both 10% scenarios lead to lower overall conservation values than the 17% connectivity-based scenario, they both largely increased the retention of the different conservation criteria in comparison with BAU (Fig. 4). The 10% quality and connectivity-based scenario was better than, or at least almost as good as, the 10% quality-based scenario for all the criteria, which implies that trade-offs in the data were not too large. The 10% quality and connectivity-based scenario was also more efficient than the 17% quality-based scenario at retaining long-range connectivity (Fig. 3) but not the other conservation criteria (Fig. 4).

DISCUSSION

We used graph-based connectivity analyses to design habitat networks for biodiversity conservation that account for uncertainties in both future land-use and climate change. To do so, we solved several computational challenges of regional-scale analyses (e.g., analyzing landscapes with several thousand of nodes [START_REF] Pelletier | Applying Circuit Theory for Corridor Expansion and Management at Regional Scales: Tiling, Pinch Points, and Omnidirectional Connectivity[END_REF][START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]) and combined recent advances in habitat-network modeling and reserve design into a single workflow. Our framework goes one step further than the previous attempts to include connectivity into spatial prioritization that were conceptual [START_REF] Lehtomäki | Methods and workflow for spatial conservation prioritization using Zonation[END_REF], static [START_REF] Magris | Integrating multiple species connectivity and habitat quality into conservation planning for coral reefs[END_REF] or single-species oriented [START_REF] Rayfield | Multipurpose habitat networks for short-range and long-range connectivity: a new method combining graph and circuit connectivity[END_REF]. Moreover, by combining fine-scaled land-use simulations with a set of regional climate projections (Coulter Riordan & Rundel 2014), we evaluated the performance and robustness of our case study habitat networks to future threats. Thus, we join the few researchers who have tackled the consequences of both land-use and climate change on biodiversity (but see, e.g., [START_REF] Jetz | Projected impacts of climate and land-use change on the global diversity of birds[END_REF][START_REF] Mazaris | Evaluating the Connectivity of a Protected Areas' Network under the Prism of Global Change: The Efficiency of the European Natura 2000 Network for Four Birds of Prey[END_REF]. We kept our framework as general as possible and applied it to a landscape that has features typical of highly fragmented urbanizing landscapes around the world. Our framework and methods can thus be applied to any other highly fragmented landscape; include other connectivity measures; include other types of conservation criteria (e.g. ecosystem services, [START_REF] Jantz | Carbon stock corridors to mitigate climate change and promote biodiversity in the tropics[END_REF]; or be derived with other parameterizations (e.g. prioritization weighting schemes). Some of our results inform land planning and biodiversity conservation under climate change. First, we showed that managing landscapes specifically for connectivity and habitat quality remains a potential lever for conservation action because it slows the breakdown of sparse habitat networks subjected to ongoing land-use change. Accounting for diverse connectivity facets in spatial prioritization strongly modified conservation priorities [START_REF] Crouzeilles | The importance of using sustainable use protected areas for functional connectivity[END_REF]) and connectivity conservation scenarios were especially effective at protecting long-range connectivity. Protecting the top 10% of quality and connectivity-based priorities was more efficient than protecting the top 10% of quality-based priorities for the retention of all conservation criteria and almost as efficient as protecting the top 17% of quality and connectivity-based priorities for the maintenance of long-range connectivity.

Nonetheless, little of the remaining forest can be lost without some loss of conservation value, and short-range connectivity (here related to the metapopulation carrying capacity) was highly sensitive to further habitat loss especially for long-distance dispersers [START_REF] Mazaris | Evaluating the Connectivity of a Protected Areas' Network under the Prism of Global Change: The Efficiency of the European Natura 2000 Network for Four Birds of Prey[END_REF]) due to the loss of large stepping-stone habitat patches [START_REF] Saura | Stepping stones are crucial for species' long-distance dispersal and range expansion through habitat networks[END_REF]). This will be true of regions, like ours, that are already heavily modified by human activities, where natural habitat area is well below thresholds at which negative impacts of fragmentation become critical [START_REF] Andren | Effects of habitat fragmentation on birds and mammals in landscapes with different proportions of suitable habitat: a review[END_REF][START_REF] Saura | Network analysis to assess landscape connectivity trends: application to European forests (1990-2000)[END_REF]. Although prioritizing a network for connectivity and habitat quality can lead to efficient outcomes for biodiversity conservation it cannot compensate for further habitat loss [START_REF] Minor | A multiscale network analysis of protected-area connecitivty for mammals in the United States[END_REF], which would require a reduction of land-conversion rates and the restoration of habitat optimally contributing to short-and long-range connectivity [START_REF] Zetterberg | Making graph theory operational for landscape ecological assessments, planning, and design[END_REF]. To identify areas where connectivity can be improved, such as the implementation of corridors to reduce dispersal bottlenecks, our framework could be tailored with criteria related to the importance of links for both short-and long-range connectivity.

Second, in our case study, spatial conservation prioritizations made without climate as a criterion were robust to future climate change out to 2050, especially the quality and connectivity-based prioritization. This indicates the uncertainties related to climate change (socioeconomic scenarios, circulation models, species distribution models) did not greatly affect the prioritization of habitat patches for conservation within this habitat network (but see [START_REF] Mumby | Reserve design for uncertain responses of coral reefs to climate change[END_REF][START_REF] Kujala | Conservation Planning with Uncertain Climate Change Projections[END_REF]). This event is likely in ecotone regions between temperate and boreal ecosystems, where species richness is expected to be driven by landscape heterogeneity more than by climate [START_REF] Kerr | The impact of climate change on mammal diversity in Canada[END_REF] and where climate is expected to change rather homogeneously [START_REF] Berteaux | The CC-Bio Project: studying the effects of climate change on Quebec biodiversity[END_REF]. In low-contrast temperate regions, conservation priorities are expected to be little modified by the interactions between land-use and climate change [START_REF] Jetz | Projected impacts of climate and land-use change on the global diversity of birds[END_REF][START_REF] Mantyka-Pringle | Climate change modifies risk of global biodiversity loss due to land-cover change[END_REF]. It would be valuable to apply our framework to other more climatically heterogeneous or extreme regions over longer periods to examine how climate change would alter connectivity-based prioritization. Also, methods to combine current connectivity with changing connectivity through time are emerging and should make it possible to optimize the prioritization and loss process [START_REF] Hodgson | How to manipulate landscapes to improve the potential for rane expansion[END_REF].

Third, we included multiscale and species-specific parameters for dispersal ability and habitat requirements, thus rooting our approach in metapopulation theory [START_REF] Saura | Network analysis to assess landscape connectivity trends: application to European forests (1990-2000)[END_REF]. We believe this is necessary to move from structural to functional connectivity assessments, which cannot be properly done with coarsely defined habitat maps and dispersal distances.

Future research on functional connectivity and connectivity dynamics needs to better connect network theory with field data and individual-based stochastic movement simulators to update resistance surfaces and dispersal kernels, validate dispersal pathways, and quantify the spatiotemporal connectivity among patches through time [START_REF] Bocedi | RangeShifter: a platform for modelling spatial eco-evolutionary dynamics and species' responses to environmental changes[END_REF]. For example, movement data associated with genetic analyses have helped improve connectivity estimates and so validated the ongoing northward range expansion of the white-footed mouse in our region (e.g. [START_REF] Marrotte | Landscape resistance and habitat combine to provide an optimal model of genetic structure and connectivity at the range margin of a small mammal[END_REF].

We identified 3 key recommendations to guide regional land planning. First, it is crucial to combine multiple criteria for both connectivity and habitat quality into a prioritization analysis for multiple species at a species-specific level (but see, [START_REF] Magris | Integrating multiple species connectivity and habitat quality into conservation planning for coral reefs[END_REF]. The relative contribution of individual patches to the regional connectivity depends on multiple factors, and only by accounting for these divergent conservation criteria can a good balance between the current and future habitat and connectivity requirements of all species within the network be ensured [START_REF] Meller | Ensemble distribution models in conservation prioritization: from consensus predictions to consensus reserve networks[END_REF]. We found that conserving connectivity and habitat quality for multiple species requires substantial increases in total habitat area protected and efficient prioritization of the habitat to be protected based on multiple criteria.

Second, we recommend the use of realistic scenario-based projections of climate and landuse change so as to future-proof the structure of the habitat network being prioritized. Given the investment required to acquire, restore, and protect habitat within a network, it is important for the network to be robust to as many likely scenarios of change as possible [START_REF] Titeux | Biodiversity scenarios neglect future land-use changes[END_REF]. Although experiments at the landscape scale are impractical [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF], scenario simulations can be used to capture the range of impacts of future fragmentation on biodiversity and explore how regional network planning will interact with regional development to modify connectivity. Landscape simulations go beyond basic graphbased simulations because land-use change does not lead to random patch deletion but rather to spatially structured changes that affect the habitat network including the matrix surrounding the protected areas.

Third, there is value in working with stakeholders from the start [START_REF] Beier | Toward best practices for developing regional connectivity maps[END_REF][START_REF] Brodie | Connecting science policy and implementation for landscape scale habitat connecitivty[END_REF]). In the case of Montreal, stakeholders have been engaged at many levels (e.g. governments, environmental nongovernmental organizations, agricultural unions, and private landowners) [START_REF] Mitchell | The Montérégie Connection Project : Linking landscapes, biodiversity, and ecosystem services to improve decision-making[END_REF] during and after the modeling process. This engagement helped bridge a disparate array of conservation initiatives (i.e. corridor projects with plans for the city greenbelt) into a regional plan for connectivity conservation. The ongoing implementation of the network has been fostered by an open process allowing for the codevelopment of the vision and goals of the ecological network. Additional criteria have since been added (e.g. ecosystem services) to broaden the value and appeal of the network to functional dimensions of biodiversity, which are vital to human well-being in the region.

We showed that accounting for connectivity in spatial prioritization strongly modifies conservation priorities. Through our framework, we applied multiple connectivity criteria to meet the habitat and dispersal needs of many species facing uncertain future land-use and climate change. The quality and connectivity-based conservation priorities we identified for periurban Montreal were area efficient and provided favorable trade-offs among multiple conservation criteria. Our framework is general and well suited to the design of ecological networks and green infrastructure for biodiversity and ecosystem services in other regions, in particular regions where habitat connectivity is critically low. These are potential criteria retention levels from Zonation without applying land-use change or recalculating network metrics. The retention of all criteria is displayed, even the climate criterion that is not targeted by this prioritization scheme. Stacked columns highlight the relative potential conservation value of criteria for our business-as-usual (BAU) (the 1.2% currently protected) and conservation scenarios (10% and 17% of area protected) and for all habitat protected (30%). The overall mean potential conservation value of each scenario is displayed in grey. The stacked column for *Quality (10%) and *Quality (17%) correspond to the quality-based prioritization scheme (Supporting Information) and are displayed here to emphasize differences between quality-and quality and connectivity-based prioritizations.

Relative change in equivalent connectivity after systematic node removal (dEC) quantifies node importance in terms of the total amount and quality of habitat species can reach (Table 2). 'Climate' is an average of the climatic suitability layers across species, climate projections and horizons (2025,2050).

  protected areas covered the 1.2% of terrestrial St. Lawrence Lowlands actually protected in 2000. Four conservation scenarios were developed within a factorial design to test the quantity of protected areas (expansion of current protected areas based on the top-ranked spatial conservation priorities to reach 10% or 17% of the area) and criteria used to define them (prioritizations based on either quality [Fig. 2a] or quality and connectivity criteria [Fig. 2b]). Climate criteria were not included in the conservation scenarios because their inclusion did not have a strong enough effect on the configuration of priority areas. We chose 17% of protected land to follow target 11 of the Aïchi biodiversity targets (Convention on Biological Diversity 2010) and 10% of protected land arbitrarily as an intermediate level of protection.

  Depending on the species, habitat of focal species covered 4-25% of the St. Lawrence Lowlands area. The habitat networks had 113-16,590 patches (mean area 35-5,989 ha), and 8 (Pileated woodpecker) to 124 (Wood frog) km of links were required for individuals to traverse the 150 km separating the Laurentian and Appalachian Mountains (Supporting Information).

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Comparison of spatial conservation prioritizations based on (a) habitat quality, (b) habitat quality and connectivity, (c) habitat quality and climate suitability, and (d) habitat quality, connectivity and climate suitability for 14 focal species (5, 10 and 17% areas with highest priorities).

Figure 3 .

 3 Figure 3. Changes in (a, b) short-range (equivalent connectivity in St. Lawrence Lowlands) and (c, d) long-range (mean cumulative resistance of shortest paths between mountain ranges) connectivity for the (a, c) business-as-usual (BAU) scenario over time and for (b, d) all conservation scenarios in 2050 (open symbols, mammals; closed symbols, birds; line symbols [e.g., +], amphibians; blue, low dispersal ability; green, yellow, orange; moderate dispersal ability; red, high dispersal ability; black dots, average effect across species). Scenario comparisons in 2050 show the effects of protecting habitat in the context of ongoing land use: BAU (1.2% of area currently protected) versus 10% or 17% of area protected based on habitat quality (10Q and 17Q) or habitat quality and connectivity (10Q&C and 17Q&C).

Figure 4 .

 4 Figure 4. The mean percentage of conservation criteria (mean among the fourteen focal species) retained by Zonation in the quality and connectivity-based conservation prioritization scheme as a function of the percentage of the St. Lawrence Lowlands protected.

Acknowledgments

This research was funded by the OURANOS ÉCOBIOCC project (project 554014) with a grant to A.G., E. Bennett, J. Cardille, and M. Lechowicz. C.H.A. was supported by a Marie Curie International Outgoing Fellowship under the 7th European Community Framework Program (DYVERSE project, no. 272284). B.R. was supported by an NSERC Postdoctoral Fellowship. A.G. was supported by a Canada Research Chair in Biodiversity Science, NSERC Discovery grant, and the Liber Ero Chair in Conservation Biology. Computations

were performed using the Amazon EC2 cloud computing platform. We thank A. Millon, D.

Renard, and A. Doxa for comments on earlier drafts and D. Pelletier for performing circuitbased connectivity analyses. We also thank A. Dabrowski, K. Dancose, L. Meulien, V.

Supporting Information

Supplementary figures and methods (Appendix S1) are available online. The authors are solely responsible for the content and functionality of these materials. Queries (other than absence of the material) should be directed to the corresponding author.

Literature Cited

Table2. Summary of habitat quality and connectivity metrics used as input criteria in the prioritization schemes.


Wolfgang Cramer

Recognizing the work of journal editors and reviewers

Recognizing the work of journal editors and reviewers

Wolfgang Cramer & James Ford Regional Environmental Change now enters its 17 th volume, having published 1161 papers in 75 issues through its existence, 209 of them in 2016 alone. As editors, we are quite proud about this. We also know just how much work this represents, first of all by the authors presenting their work to us, and then through the huge time investment and intellectual contribution made by subject editors, guest editors and reviewers. Much of this work is carried out in anonymity, and so authors have little opportunity to thank for the invaluable support of their work they are receiving through this process. But we know this work is appreciated! For the papers submitted during the period of Nov 1, 2015 to Oct 31, 2016, in total 387 reviewers provided advice, often in the form of detailed and constructive statements including specific remarks for parts of the manuscript. On behalf of the authors, we as editors wish to thank these colleagues, and also all other reviewers, for their hugely important contribution. From next year onwards we plan to give this work some recognition by listing the full list of reviewers in the last issue of any year.

We are pleased to be able to work with a strong group of handling editors -also they contribute large amounts of their time to the common objective of producing publications of high academic quality and interdisciplinary standards. Three of them deserve to be mentioned specifically, since they are now stepping down from work for our journal, after many years of unfaltering support: Elena Bennett, Helmut Geist, and Will Steffen have each been handling many papers, engaged in discussions with reviewers, authors and chief editors, and generally sorted out numerous small and large problems before papers were published. We want to use this space to point out just how grateful we are for their support! Now, with the first issue of 2017, we have decided to re-organize the journal's crucial support structure, its Editorial Board. In the past, the board consisted of a large group of colleagues (64 at the time of writing) who had all supported us by a great number of reviews and other activities. From now on, we plan to make their contributions visible by listing all reviewers at the end of a given annual volume. The Editorial Board, from now, on, will be composed of only the active handling editors, a group that is growing in number and also in diversity across disciplines and countries of residence. It currently includes Robbert Biesbroek, Virginia Burkett, Erika Coppola, Nicolas Dendoncker, Xiangzheng Deng, Sarah Gergel, Helmut Haberl, Juan Ignacio Lopez Moreno, Shuaib Lwasa, Ülo Mander, Marc J. Metzger, Jamie Pittock, Christopher Reyer, Diana Sietz, and Erica Smithwick -but several more colleagues will join us in the near future. So we thank the old Editorial Board for its reviewing activities (which we hope will continue), and we welcome our new Editorial Board with whom we plan to engage also in terms of discussing more generic issues of journal development.

Finally, the end of 2016 also has marked the end of the involvement of the journal's most important person for daily management and continuity. Our Managing Editor Ms. Gabriele Götz has been the key person for virtually all communications between authors, editors and reviewers, and she has set up and managed our extensive archiving system. But most of all, Gabriele Götz has developed the "art of reminding people". We are writing this with no ironic undertones: it is clear that, these days, most academics work under high pressure in a complex multi-tasking environment. Reminding them of earlier commitments or invitations, such as to review a paper, or for an editor, to invite more reviewers, is not just a technical issue, it is also a human interaction that can be managed with a smile and a personalised way of considering every person's specific situation. Gabriele Götz has carried this task over fifteen years, interacting with thousands of individuals, and her share in the success of Regional Environmental Change can simply not be overestimated. Gabi, all editors will miss the exchanges with you, and we hope many authors and reviewers have also understood how helpful your work has been for them. At this point, the editors also wish to thank the Potsdam Institute for Climate Impact Research to have hosted the Editorial Office, and provided the staff time for Gabriele Götz over this entire period.

Concluding, we want to state briefly our objectives for 2017 and the forthcoming years. First, we will carry on with our Editorial Board and our Guest Editors in producing what we hope are excellent scientific contributions to the inter-and transdisciplinary theme of regional environmental change. We will do this by expanding the journal's Editorial Board and by continuing to consider strong proposals for Special Issues edited by Guest Editors. While we see a continuing role of 12 page full articles as before, we also want to encourage "rapid communications" in order to inject particular ideas or findings into the scientific literature. On a technical level, we hope to return to smooth operations by establishing a new editorial office with a new managing editor in the near future -the journal's website and front matter will report progress about this as soon as this has been achieved.

It is an inevitable aspect of success for a journal that many submitting authors may not see their papers eventually be published, this general situation will not change. During the years 2013-2015, the rejection rate of the journal has fluctuated between 60 and 70%. We wish to point out, however, that more than half of the papers which are not published in Regional Environmental Change are invited for transfer elsewhere, often accompanied by constructive reviews. Many rejections could have been avoided by careful study of our journal's goal statement. We continue to strive to fair handling of all submissions, independent of their origin or the experience of the submitting authors.