Jesús Díaz-Sanz Ab

Samuel Robert

Catherine Keller

Jesús Díaz-Sanz

Parameters influencing run-off on vegetated urban soils: A case study in Marseilles, France

Keywords: Infiltration, Compaction, Run-off, Urban soils, Technosols, Tile industry

In urban areas, episodes of extreme precipitation may generate run-off floods, but urban soils are seldom considered in spatial planning documents. Our aim was to assess the infiltration of vegetated urban soils and to identify the soil parameters affecting infiltration in order to explore the potential of these soils for run-off mitigation under a Mediterranean climate with extreme precipitation episodes. We compared vegetated soils in the "urban" zone of the Marseilles' 16 th district (France) either lying on former clay quarries or tile factories (T) or not (NT) and with different land-uses. We used a simplified method to measure field saturated hydraulic conductivity (Kfs), measured penetration resistance in the first 0.15 m (Qd) and analysed physico-chemical properties in the first horizon. Levels of runoff production were estimated according to Kfs and an adaptation of USDA Hydrologic Soil Groups. All soils were Technosols or Anthrosols regardless of their prior history. We found that land-use history influenced Kfs, with T soils and soils with current vegetated land-use younger than 13 years having lower Kfs than NT soils and soils with land-use older than 13 years. Soil organic matter content influenced Kfs positively when soils were young, while vegetation cover and penetration resistance (most of the soils had Qd > 2 MPa) had not a clear effect on Kfs. Overall, combining Kfs and soil depth, 14.3 % of the T soils and 84.6 % of the NT soils had low levels of run-off production. It is therefore recommended that 1) Technosols or Anthrosols should be characterised for their unique physicochemical and physical properties but also for their land-use history, and 2) the infiltration of these soils should be considered in spatial planning documents.

Introduction

Extreme precipitation episodes may affect urban areas. The intensity of these episodes and poor infiltration in soils may generate run-off floods [START_REF] Bhaduri | Long-term hydrologic impact of urbanization: A tale of two models[END_REF][START_REF] Chalvet | Marseille, une ville méditerranéenne entre pénurie et inondation[END_REF][START_REF] Kundzewicz | Flood risk and climate change: global and regional perspectives[END_REF]. These floods produce substantial damage to the urban network, increase traffic jams and accidents, and are a source of pollution [START_REF] Davis | Loading estimates of lead, copper, cadmium, and zinc in urban runoff from specific sources[END_REF][START_REF] Koetse | The impact of climate change and weather on transport: An overview of empirical findings[END_REF][START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF]. Several factors increase the risk of run-off flooding: (1) a higher intensity of precipitations [START_REF] Trenberth | Changes in precipitation with climate change[END_REF] and ( 2) urbanization [START_REF] Kundzewicz | Flood risk and climate change: global and regional perspectives[END_REF]. Urbanization jeopardizes the territory, reducing permeable surfaces such as those found in French Mediterranean coastal zones [START_REF] Robert | A framework to analyse urban sprawl in the French Mediterranean coastal zone[END_REF].

Urbanization also affects infiltration in soils since it contributes to compaction [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF][START_REF] Shuster | Residential demolition and its impact on vacant lot hydrology: Implications for the management of stormwater and sewer system overflows[END_REF].

One of the goals of planning documents is to tackle flood risks in order to reduce both hazards and vulnerability. Water management strategies encourage the collection of stormwater to alleviate the hazard, either by adapting the sewer network or using sustainable drainage systems [START_REF] Perales-Momparler | A regenerative urban stormwater management methodology: The journey of a Mediterranean city[END_REF]. Spatial planning may also limit urbanization, but local planning practices can be very uneven, as shown in Southern France [START_REF] Prévost | Local spatial planning practices in four French Mediterranean coastal territories under pressure[END_REF]. Generally, urban planning documents do not explicitly consider urban soils as a means to minimize run-off. In the case of France, local spatial planning documents -Plan Stormwater infiltration relies generally on the characteristics of (1) soil and (2) ground surface. The major parameters are: (1) soil texture, structure, humidity, depth, soil organic carbon (SOC), compaction and existence of impervious layers [START_REF] Elrick | Hydraulic conductivity measurements in the unsaturated zone using improved well analyses[END_REF][START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Mayor | Factors and interactions controlling infiltration, runoff, and soil loss at the microscale in a patchy Mediterranean semiarid landscape[END_REF]; and (2) slope, vegetal cover type and density, presence of rock fragments, artifacts, cracks, biological activity or crusts on the ground surface [START_REF] Mayor | Factors and interactions controlling infiltration, runoff, and soil loss at the microscale in a patchy Mediterranean semiarid landscape[END_REF][START_REF] Mavimbela | Characterizing infiltration and internal drainage of South African dryland soils[END_REF][START_REF] Morbidelli | Role of slope on infiltration: A review[END_REF]. Urban areas have a mosaic of relict natural soils and other soils, strongly impacted by anthrosolization (human activities). The removal of topsoil, backfilling with fine and coarse materials, earthworks and ploughing lead to Anthrosols (AT) and Technosols (TC). According to the IUSS Working Group WRB (2015), AT are soils generated by intensive agricultural landuse and backfilling with fine materials, while TC contain coarse materials produced by humans (artefacts and technical material). Anthrosolization thus generates different hydrological properties in contrast to natural soils [START_REF] Séré | Predictability of the evolution of the soil structure using water flow modeling for a constructed Technosol[END_REF].

The presence of coarse materials in the soil matrix of TC diminishes the area available for infiltration [START_REF] Yilmaz | Physical properties of structural soils containing waste materials to achieve urban greening[END_REF]. In addition, compaction in AT and TC limits infiltration since it reduces the bulk density and macroporosity [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Shuster | Residential demolition and its impact on vacant lot hydrology: Implications for the management of stormwater and sewer system overflows[END_REF]. The impact of compaction in these soils is greater than in natural soils because they already present a limited porosity [START_REF] Grabosky | Plant available moisture in stone-soil media for use under pavement while allowing urban tree root growth[END_REF]. Furthermore, the backfill materials used may induce specific processes such as precipitation of carbonates that can narrow the pore size [START_REF] Wissmeier | Effect of mineral reactions on the hydraulic properties of unsaturated soils: Model development and application[END_REF][START_REF] Yilmaz | Influence of Carbonation on the Microstructure and Hydraulic Properties of a Basic Oxygen Furnace Slag[END_REF] and subsequently the water storage, exposing vegetation to water stress [START_REF] Urban | Two Different Approaches to Improve Growing Conditions for Trees Comparing Silva Cells and Structural Soil[END_REF][START_REF] Huot | Lysimeter monitoring as assessment of the potential for revegetation to manage former iron industry settling ponds[END_REF].

Although the nature and composition of AT and TC can limit infiltration, some of them can significantly reduce run-off [START_REF] Xiao | Performance of engineered soil and trees in a parking lot bioswale[END_REF]. Coarse materials can create fissures at the interface with fine backfill materials [START_REF] Cousin | Influence of rock fragments on the water retention and water percolation in a calcareous soil[END_REF][START_REF] Wang | Effects of Urbanization, Soil Property and Vegetation Configuration on Soil Infiltration of Urban Forest in Changchun, Northeast China[END_REF]. Similarly, vegetation cover and organic matter can enhance porosity [START_REF] Layman | Below ground matters: Urban soil rehabilitation increases tree canopy and speeds establishment[END_REF][START_REF] Wiesner | Spatial and temporal variability of urban soil water dynamics observed by a soil monitoring network[END_REF].

The study of hydraulic properties on AT and TC can be done with infiltration techniques originally developed for natural soils. This is the case for methods based on the Beerkan estimation of soil transfer parameters (BEST), which were used to calculate hydraulic properties accurately on TC (e.g. the field-saturated soil hydraulic conductivity (Kfs), the soil water retention curve or the saturated hydraulic curve) [START_REF] Yilmaz | Assessment of hydraulic properties of technosols using Beerkan and multiple tension disc infiltration methods[END_REF]. Nevertheless, BEST methods require the determination of the bulk density, the soil water content prior to the experiment and the particle size distribution [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through Infiltration Experiments-BEST[END_REF][START_REF] Yilmaz | Assessment of hydraulic properties of technosols using Beerkan and multiple tension disc infiltration methods[END_REF]. The numerous requisites of the BEST methods made difficult to study hydraulic properties in every field condition. So, simplified approaches were developed to estimate some hydraulic properties that were mainly limited to Kfs. Some examples of simplified approaches are the One-Ponded Height Technique [START_REF] Elrick | Hydraulic conductivity measurements in the unsaturated zone using improved well analyses[END_REF], the Simplified Beerkan Infiltration test (SBI) [START_REF] Bagarello | Estimating fieldsaturated soil hydraulic conductivity by a simplified Beerkan infiltration experiment[END_REF] or the Wu1 and Wu2 methods [START_REF] Wu | Measuring saturated hydraulic conductivity using a generalized solution for single-ring infiltrometers[END_REF][START_REF] Xu | Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of BEST and Wu methods[END_REF], which have not been applied to AT or TC soils so far. On natural vegetated soils of different textures, Wu1 had similar Kfs results to BEST_slope [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through Infiltration Experiments-BEST[END_REF] and Wu2 to BEST_intercept (Yilmaz et al., 2010), but both Wu methods permitted to calculate Kfs in more soils independently of the texture [START_REF] Xu | Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of BEST and Wu methods[END_REF].

In general, the contribution of AT and TC to stormwater infiltration and thus to the mitigation of run-off in urban areas has not been fully assessed yet [START_REF] Grabosky | Plant available moisture in stone-soil media for use under pavement while allowing urban tree root growth[END_REF][START_REF] Fletcher | Understanding, management and modelling of urban hydrology and its consequences for receiving waters: A state of the art[END_REF]: the type of soil is often not detailed, preventing comparison and/or generalization of the results, and Mediterranean urban areas, prone to intense episodes of precipitation, have been little studied. Moreover, the impact of the land-use history needs to be clarified, even if some studies showed more infiltration on urban soils with older age of land-use [START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF][START_REF] Shuster | Residential demolition and its impact on vacant lot hydrology: Implications for the management of stormwater and sewer system overflows[END_REF].

For these reasons, the objectives of this study were (1) to assess the potential of vegetated urban soils to mitigate run-off in an area prone to intense episodes of precipitation via infiltration studies, (2) to identify the soil parameters that affect infiltration, and (3) to assess the influence of the land-use history on infiltration and, more specifically, the industrial history (tile industry). Field measurements and sampling were performed in Marseilles (France), in plots covering a large range of known current and past land-uses.

Materials and methods

Study area

Marseilles (France) is a coastal city of the Mediterranean Sea, with a little more than 860,000 inhabitants (Insee, 2019). Marseilles is a landmark for the French economy due to its port and industrial activities [START_REF] Prévost | Local spatial planning practices in four French Mediterranean coastal territories under pressure[END_REF] The study area (Fig. 1) is in the "urban" zone of the 16 th district, as defined and delimited in the urban planning document, Plan Local d'Urbanisme (PLU) [START_REF] Mamp | Plan Local d'Urbanisme de Marseille -Rapport de présentation[END_REF]. The 16 th district marks the northern city limit, with a transition towards the surrounding natural areas. The urban fabric is strongly influenced by the commercial harbour and the associated industrial and logistical zones [START_REF] Robert | Entre étalement et densification : une approche fine de l'urbanisation littorale sur la Côte bleue, Provence[END_REF]. The "urban" zone extends over 7.2 km 2 and includes 1.0 km 2 of semi-natural and green zones, 2.9 km 2 of economic zones, 1.1 km 2 of transportation infrastructure, and 2.1 km 2 of low density residential areas (data drawn from the "Occupation des sols à grande échelle 2011 -Côte Bleue" Database, developed and presented in [START_REF] Robert | Entre étalement et densification : une approche fine de l'urbanisation littorale sur la Côte bleue, Provence[END_REF]. This zone is currently undergoing a process of urban renewal that jeopardizes unsealed soils [START_REF] Robert | Entre étalement et densification : une approche fine de l'urbanisation littorale sur la Côte bleue, Provence[END_REF]. A particularity of the study area is that it has an industrial past dominated by the tile industry. A sampling strategy was therefore designed to allow for comparison between the area impacted by this industry (located in a clay basin) and the rest of the "urban" area, which is assumed not to be impacted by this activity.

The area impacted by the tile industry covers circa 16.7 % of the "urban" zone, a percentage drawn from photointerpretation of the 1957 aerial image [START_REF] Ign | C3245-0461_1957_CDP1202[END_REF]. It includes former clay quarries for extraction, and tile factories. This area underwent intense earthworks and was backfilled with construction debris (bricks, tiles, gravel).

Spatial data base construction

A geographical database comprising aerial images, topographic maps, digital elevation models, geological data, land use maps, roads and railways, buildings, land planning zones and other data relevant to the study was organised into a GIS, using ArcGIS (Fig. 2

Soil data acquisition

Authorization was obtained to sample vegetated urban soils from June 2017 to June 2018: (1) 14 plots lying on former clay quarries and tile factories (T); (2) 13

plots not lying on former clay quarries or tile factories (NT); (3) one control plot in the surrounding shrubland not subjected to urbanization (Fig. 1). Given the sampling restrictions on urban soils, we performed one-shot sampling in the representative area of each plot.

In the field, we categorized soil surfaces according to the vegetation cover as trees only, herbaceous & trees, herbaceous only, and herbaceous & shrub. Then we described soil profiles from auger cores or pit observations where possible. In the description, we noted soil depth, and described and sampled the first horizon (average thickness of 0.13 m).

Penetration resistance measurement

In addition to the sampling, we measured in situ the penetration resistance (Qd)

and infiltration. Qd gives information on soil compaction [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Keller | Soil Quality Assessment for Spatial Planning in Urban and Peri-Urban Areas -Municipalities of Gardanne and Rousset (southern France)[END_REF][START_REF] Bean | Effect of Amendment Type and Incorporation Depth on Runoff from Compacted Sandy Soils[END_REF]. Qd was measured with a constant weight light dynamic cone penetrometer (Type 2, Sols Mesures, France), and the analysis was stopped if 3 equal readings were taken or if the penetrometer reached the profile depth. Then, we calculated Qd for the first 0.15 m of the soil profile since the highest compaction occurs in the first 0.15 m [START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF][START_REF] Bean | Effect of Amendment Type and Incorporation Depth on Runoff from Compacted Sandy Soils[END_REF]. To check this, Qd was also calculated for the surface horizon and for the depth profile. When there were several horizons within the 0.15 cm, Qd was weighted according to the thickness of these horizons.

Infiltration experiment

Between June 2017 and June 2018, we carried out a constant-head infiltration experiment, following the Wu2 method [START_REF] Wu | Measuring saturated hydraulic conductivity using a generalized solution for single-ring infiltrometers[END_REF][START_REF] Xu | Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of BEST and Wu methods[END_REF]. Given the restrictions on operating in urban soils, infiltration was measured once in the field, and the experiment was repeated in 11 plots to check the measurement (Supplementary Table 1). In the plots where the infiltration was repeated, the first result was used. The Wu2 method was conducted using a single cylinder with a shaped rim and an internal radius of 0.05 m. Before starting the experiment, the vegetation was cut. Then, the cylinder was driven around 0.05 m into the soil, recording the exact cylinder depth and preserving the surface to avoid lateral flow.

Water was poured onto the soil inside the cylinder until a water head of 0.10 m.

From this moment on, the time started to be measured. When the water head reached 0.09 m, water was added until reaching 0.10 m again, and the added volume and the time were recorded. This process was repeated for 2 to 3 hours until a steady state was reached. At the end of the experiment, the volumes were plotted as a function of time to obtain infiltration curves in order to identify the steady infiltration of water, which was called infiltration rate. At this step, we discarded one plot out of 29 because of infinite infiltration. The cumulative height of infiltrated water was also plotted with time to produce the curve required by the Wu2 method.

The Wu2 method assumes soil got saturated during the infiltration experiment.

According to this assumption, the cumulative infiltration curve has an initial transient part and a final steady part. Those are determined since the former has a variable slope, while the latter a constant slope. The field-saturated soil hydraulic conductivity (Kfs) is then estimated from (Eq.1):

Equation 1. Equation to estimate the field-saturated soil hydraulic conductivity (Kfs) equation according to the Wu2 method [START_REF] Xu | Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of BEST and Wu methods[END_REF].

where i is the slope of the steady part of the cumulative infiltration curve, a an adimensional parameter equal to 0.9084 and f a correction factor. The f factor can be estimated from (Eq. 2):

Equation 2. Correction factor [START_REF] Wu | Measuring saturated hydraulic conductivity using a generalized solution for single-ring infiltrometers[END_REF].

where H is the water head, α a soil texture evaluation and G* the adjustment parameter of the infiltrometer characteristics. G * can be determined by Eq. 3, where d is the cylinder depth and r the radium of the cylinder:

Equation 3. Calculation of the adjustment parameter G* according to the characteristics of the ring infiltrometer [START_REF] Wu | Measuring saturated hydraulic conductivity using a generalized solution for single-ring infiltrometers[END_REF].

Eq. 2 and 3 allow estimating Kfs in Eq. 1, but Eq.2 requires a field evaluation of the soil texture to choose α [START_REF] Wu | Measuring saturated hydraulic conductivity using a generalized solution for single-ring infiltrometers[END_REF]. The Wu parameter α takes 3 different values depending on the texture of the studied soil (Table 1). The texture class of the first horizon was assigned according to the particle-size distribution of the USDA triangle (Soil Science Division Staff, 2017). The selection of an erroneous α value does not compromise Kfs estimation [START_REF] Bagarello | Estimating fieldsaturated soil hydraulic conductivity by a simplified Beerkan infiltration experiment[END_REF]. In fact, an improper choice produces an acceptable error [START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF].

Estimation of levels of run-off production

Spatial planning requires an assessment of run-off production in urban soils to manage the risk of run-off floods. The levels of run-off production can be estimated through the combination of soil depth and Kfs to assist in the management of run-off flood risk in urban areas. Soil depth is related to the potential stormwater storage while Kfs indicates the stormwater infiltration. We assessed the level of run-off production based on an adaptation of the Hydrologic Soil Groups (HSG) (USDA-NCRS, 2009) (Table 2). Our criteria were: (1) soils have an estimated low run-off when Kfs > 10 µm s -1 and soil depth ≥ 0.50 m (HSG A and B); (2) soils have an estimated high run-off when Kfs ≤ 10 µm s -1 (HSG C and D), or when depth is < 0.50 m down to the impervious layer or to a very dense material (e.g. technogenic material, bedrock) (HSG D). We did not quantify run-off production since our objective was to assess the infiltration of vegetated urban soils to mitigate run-off.

Furthermore, we did not test the influence of slope since most of the urbanized area and the sampled plots were flat or on terraces.

Physico-chemical analyses

The soil samples (first horizon) were dried at 40ºC and sieved at 2 mm and further ground when needed. Samples were analyzed for soil organic carbon (SOC) (AFNOR, 1995a), total nitrogen (N) and Olsen phosphorus (P) (AFNOR, 1995b(AFNOR, , 1995c)), and particle-size distribution [START_REF] References | Qualité des sols -Détermination de la distribution granulométrique des particules du sol -Méthode à la pipette[END_REF], as these factors may influence the infiltration experiment [START_REF] Elrick | Hydraulic conductivity measurements in the unsaturated zone using improved well analyses[END_REF][START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Mayor | Factors and interactions controlling infiltration, runoff, and soil loss at the microscale in a patchy Mediterranean semiarid landscape[END_REF].

Soil data statistical treatment

The statistical analyses of soil data were performed in R (R Core Team, 2017).

Prior to the significance test, normality was checked using the Shapiro test, QQ plots and density plots. Kfs, Qd and infiltration rate had a non-normal distribution.

We investigated the significance of mean differences in (1) infiltration rate, Kfs, Qd and soil depth to assess the studied soils, and (2) Kfs to study the influence of the age of the current land-use and the vegetation cover on Kfs. Owing to the non-normal distribution of the parameters, non-parametric multivariate statistics were used, namely the Kruskal-Wallis test and a post-hoc Dunn's test with the packages FSA, car and rcompanion [START_REF] Mangiafico | rcompanion: Functions to Support Extension Education Program Evaluation[END_REF][START_REF] Ogle | FSA: Fisheries Stock Analysis[END_REF][START_REF] Fox | An {R} Companion to Applied Regression[END_REF]. The parameters affecting Kfs were explored with the Pearson test using the packages Hmisc and Corrplot [START_REF] Harrell | Hmisc: Harrell Miscellaneous[END_REF][START_REF] Wei | R package "corrplot[END_REF].

Results

Urban soil characteristics

The studied urban soils were predominantly Technosols, followed by Anthrosols (Table 2 Supplementary material). The coarse materials were composed of construction debris, gravel for all soils and additional brick fragments in T soils.

The only natural soil was a Regosol (C) (IUSS Working Group WRB, 2015) found in the surrounding shrubland. Another Regosol was found in one NT plot, in a private garden, where it was considered as a relict soil.

Soils were thicker in urbanized areas than in the C area (Fig. 3). A majority of soils had a depth exceeding 0.50 m. Although the difference was not significant, NT soils were deeper than T soils.

SOC, N and P contents

Anthrosols and Technosols had a high content of soil organic carbon (SOC), total nitrogen (N) and available phosphorus (P). Although the T soils had less SOC, N and P than NT soils, the variability of these parameters was too high (Table 3) to detect any significant difference between the two groups. In general NT soils had more organic matter and T soils had less SOC and N than the C Regosol. Moreover, the NT Regosol contained twenty times more P, 1.7 times more SOC and 1.4 times more N than the C Regosol.

Mechanical properties

The mechanical properties of urban soils may be relevant for the infiltration of stormwater and hence for the run-off production of soils. The particle-size distribution of the groups of studied soils varied slightly, but non-significantly (Fig. 4). T soils had more clay and silt than NT and C soils while NT and C soils had the same clay content, and NT contained less silt than the C one. Nevertheless, most of soils had a loam soil texture in the first horizon (Table 2 Supplementary material).

Qd for the first 0.15 m was > 2 MPa for most of the studied soils, including 100 % of the T soils, and was significantly higher in T than in NT soils (Fig. 3).

Levels of run-off production of urban soils

The studied Anthrosols and Technosols were able to infiltrate stormwater. The infiltration was similar in the NT and T plots, but lower than in the C plot (Fig. 3).

The infiltration curves showed that T and NT soils experienced a strong variability (Fig. 5). The infiltration experiment allowed for the estimation of the field-saturated hydraulic conductivity (Kfs). The Kfs relates to the infiltration of urban soils that can be estimated when the curve of the cumulative infiltration reaches a steady state. At saturation, it is possible to apply the Wu2 linear equations to estimate Kfs (Fig. 6).

The Kfs coupled with the soil depth enabled the estimation of the levels of run-off production on urban soils to be estimated. Our results highlighted that some vegetated Anthrosols and Technosols could mitigate run-off. We found that 14.3 % of T and 84.6 % of NT soils had an estimated low run-off production (Table 3). The estimated run-off production was not modified when the infiltration experiment was replicated (Table 3 Supplementary material). Since the run-off production relied both on the Kfs and the soil depth, the latter being constant for a given plot, the only parameter that could modify the estimation of run-off production in alternate infiltration experiments was Kfs.

The results of Kfs were similar to the results of infiltration rates, showing that T and NT soils had lower Kfs than C (Fig. 3). In addition, the geometric mean of Kfs in T was two times lower than in NT soils (Table 3). Furthermore, Kfs was not correlated to Qd, the particle-size distribution, or to P content in both T and NT soils (Table 4). There was a trend between Kfs and SOC and N in both types of soils, although not significant. Kfs was only significantly correlated to N in T soils (r = 0.55 and p-value = 0.039). However, N was significantly correlated to SOC in T (r = 0.97 and p-value < 0.001) and NT soils (r=0.92 and p-value < 0.001).Also, the geometric mean of SOC and N was lower in T soils (Table 3).

On the other hand, there was a significant difference in Kfs between soils with a current < 13-year-old and > 13-year-old land-use (Fig. 7). The lowest Kfs was in soils with a < 13-year-old land-use, which also corresponded to T soils. There was no significant difference in Kfs with respect to the type of existing vegetation cover, however (Fig. 8).

Discussion

Urban soil characteristics

The studied soils were mainly Technosols, with only some Anthrosols, regardless of whether they lay on former clay quarries, tile factories (T) or not (NT).

The historical presence of the tile industry induced a lower organic matter content in T soils than in NT soils, probably because the formers were younger, with a current vegetated land-use of less than 13 years. In the literature a lower organic matter content was found to be associated to the presence of brick in the coarse materials, which was also the case in the T soils. Technosols containing brick in both fine earth and skeleton had 1.8 % SOC and 1.4 g N kg -1 [START_REF] Nehls | Contribution of bricks to urban soil properties[END_REF] and 2.1 % SOC and 1.2 g N kg -1 [START_REF] Schleuß | Variability of soils in urban and periurban areas in Northern Germany[END_REF], while Anthrosols and Technosols without brick in the coarse material had 4.7 ± 3.5 % SOC and 2.4 ± 2.7g N kg -1 [START_REF] Schleuß | Variability of soils in urban and periurban areas in Northern Germany[END_REF], and 4.7 ± 3.7 % SOC and 3.3 ± 1.9g N kg -1 [START_REF] Lorenz | Biochemical characterization of urban soil profiles from Stuttgart, Germany[END_REF]. These studies also revealed a huge variability in SOC and N contents. The variability was also high in our results but T and NT soils had more SOC than the cited studies.

One characteristic of urban soils is their depth. Soil depth controls the volume of stormwater that a soil can infiltrate. Urbanization processes, such as pavement construction, can remove surficial horizons of urban soils making them shallower than other soils [START_REF] Kida | Influence of asphalt pavement construction processes on urban soil formation in Tokyo[END_REF]. However, most of the studied Anthrosols and Technosols had the required depth to store enough stormwater and avoid run-off. In our study, one of the soils which did not satisfy this requirement was the C (control) Regosol. The reason was that this soil was on a slope, unlike the N and NT soils, and had probably experienced erosion in the past, like most Mediterranean soils. In other studies, Anthrosols and Technosols were twice as deep and were never shallower than 0.50 m [START_REF] Schleuß | Variability of soils in urban and periurban areas in Northern Germany[END_REF][START_REF] Lorenz | Biochemical characterization of urban soil profiles from Stuttgart, Germany[END_REF]. The soil depth of urban soils should therefore be locally assessed to estimate the run-off production in a given urban area and cannot be estimated from the depth of natural soils.

Given the impact of human activities on urban soils, the penetration resistance for the first 0.15 m (Qd) is a critical characteristic of these soils and a proxy of soil compaction [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Keller | Soil Quality Assessment for Spatial Planning in Urban and Peri-Urban Areas -Municipalities of Gardanne and Rousset (southern France)[END_REF][START_REF] Bean | Effect of Amendment Type and Incorporation Depth on Runoff from Compacted Sandy Soils[END_REF]. The Qd for most of the studied Anthrosols and Technosols was > 2 MPa, including the totality of T soils and the C Regosol. This threshold indicates that the soils were compacted, limiting root growth [START_REF] Materechera | Influence of root diameter on the penetration of seminal roots into compacted subsoil[END_REF]Soil Science Division Staff, 2017). The C Regosol was also found to be compacted because it had numerous rock fragments originating from the bedrock. The compaction of T soils highlighted the strong impact of the earthworks associated to the demolition of tile factories and the backfilling of clay quarries.

Vegetated urban soils mitigate run-off

The infiltration of vegetated urban soils is usually not evaluated because of their supposed lower infiltration in contrast to natural soils. The decrease in infiltration was demonstrated in a study which intentionally compacted vegetated natural soils with construction machines [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF]. Similarly, in our results Kfs was lower in Anthrosols and Technosols than in the C Regosol.

Unfortunately, it is often considered that all urban soils have the same properties and a limited infiltration. Consequently, this aspect is not taken into consideration in spatial planning. However, our results highlighted that vegetated Anthrosols and Technosols can infiltrate enough stormwater to limit run-off. The methodology considered that soils generate run-off when the infiltration and the storage capacity are overwhelmed by the amount of stormwater. Both parameters, Kfs and soil depth were respectively estimated and measured for each plot. On the contrary, other studies have simulated run-off production based on Kfs only. For example, Shuster et al. ( 2014) simulated the run-off production of Technosols that had experienced demolition. They reported that soils with low field-unsaturated hydraulic conductivity (Kfun) produced as much as twenty times the run-off of soils with high Kfun (75,000 l y -1 versus 1130 to 3800 l y -1 ). The simulations of [START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF], which employed several intensities of rainfall, also found that Anthrosols and Technosols with a reduced infiltration produced more run-off. None of these studies considered the soil depth in the simulation of run-off. In our results, 66.7 % of Anthrosols and Technosols had an elevated Kfs (> 10 µm s -1 ) according to Soil Science Division Staff (2017) but only 48.1 % of these soils had low run-off (Table 3). In addition, the proportion of soils with low run-off in T soils was five times less than in NT soils. For this reason, the parameters that influence the run-off production need to be investigated more in detail.

Soil parameters affecting field-saturated hydraulic conductivity

Run-off production depends mainly on two independent factors, the soil depth and Kfs. The relevance of soil depth is that it is associated to the reservoir of stormwater that Anthrosols and Technosols can store, as previously discussed. Kfs presented a large range of values in both T and NT soils, showing that the infiltration experiment employed was useful to discriminate urban soils. Kfs variation is produced by several factors. Kfs may be sensitive to the materials contained in Technosols. The Kfs geometric mean of T soils was four times as high as that of Technosols constructed with a mixture of paper sludge [START_REF] Séré | Predictability of the evolution of the soil structure using water flow modeling for a constructed Technosol[END_REF] and two times that of Technosols composed of waste from the iron industry [START_REF] Huot | Lysimeter monitoring as assessment of the potential for revegetation to manage former iron industry settling ponds[END_REF]. In other cases, such as in Technosols containing rubble, Kfs was highly variable, ranging from 5.6 to 66.7 µm s -1 ) [START_REF] Scalenghe | The First Forty Years of a Technosol[END_REF][START_REF] Shuster | Residential demolition and its impact on vacant lot hydrology: Implications for the management of stormwater and sewer system overflows[END_REF][START_REF] Yilmaz | Physical properties of structural soils containing waste materials to achieve urban greening[END_REF].

The method employed to measure Kfs also influences the value. In [START_REF] Yilmaz | Assessment of hydraulic properties of technosols using Beerkan and multiple tension disc infiltration methods[END_REF], Beerkan estimation of soil transfer parameters (BEST) and the method with multiple tension disc infiltration, applied to a constructed Technosol composed of rubble and green waste, gave respectively Kfs values of 8.8 and 26.3 µm s -1 .

Furthermore, seasonality may influence Kfs, which appeared to be reduced when the soil was more humid [START_REF] Wiesner | Spatial and temporal variability of urban soil water dynamics observed by a soil monitoring network[END_REF]. Nevertheless, we did not observe such variation between the different campaigns because the soil was always dry according to the soil profile description.

As seen previously, compaction is a characteristic of urban soils. Compaction tends to decrease stormwater infiltration in Anthrosols and Technosols. Compaction can influence the mechanical properties of urban soils, increasing the bulk density. The higher the bulk density, the shorter the macropores and the smaller the pore size [START_REF] Gregory | Effect of urban soil compaction on infiltration rate[END_REF][START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF]. Although most of the Anthrosols and Technosols that we studied were compacted, some of them had elevated Kfs values and a low run-off production and Qd was not significantly correlated to Kfs.

density in Technosols with a construction past. Considering our results, compaction may not be a general proxy of Kfs in urban soils. We did not study porosity, but the macropore density could be responsible for the elevated Kfs values. In a study of agricultural soils, Kfs was correlated to elongated transmission pores (50-500 µm)

and fissures (> 500 µm) [START_REF] Pagliai | Soil structure and the effect of management practices[END_REF]. Provided compaction reduces elongated transmission pores in agricultural soils [START_REF] Alaoui | Does soil compaction increase floods? A review[END_REF][START_REF] Pagliai | Soil structure and the effect of management practices[END_REF], fissures may enhance Kfs in compacted urban soils. The presence of coarse material in Technosols (e.g. gravel) could create fissures at the interface with the fine material [START_REF] Cousin | Influence of rock fragments on the water retention and water percolation in a calcareous soil[END_REF][START_REF] Yilmaz | Physical properties of structural soils containing waste materials to achieve urban greening[END_REF], especially during shrinking and swelling in clay-rich soils such as the T soils. In addition, the porosity could be provided by porous coarse materials such as brick fragments [START_REF] Nehls | Contribution of bricks to urban soil properties[END_REF].

This aspect will require in-depth study in the future.

We also found that the organic matter content was only linked to Kfs in recent vegetated Anthrosols and Technosols. Kfs was significantly correlated to the N content in recent vegetated soils (<13 years), although N and SOC compose the organic matter and were significantly correlated in all soils. Recent vegetated soils (<13 years) had significantly lower Kfs and had less SOC and N. Accordingly, organic matter may be a discriminating parameter at the beginning of soil development only. At this stage, it enhances vegetation installation, root proliferation [START_REF] De Lucia | Effect of different rates of composted organic amendment on urban soil properties, growth and nutrient status of three Mediterranean native hedge species[END_REF][START_REF] Layman | Below ground matters: Urban soil rehabilitation increases tree canopy and speeds establishment[END_REF] and the formation of aggregates. These parameters increase the macropore density [START_REF] Vidal-Beaudet | Mechanical behaviour of silty clay loam/peat mixtures: cyclic compression-release tests and effects of initial water content[END_REF][START_REF] Chen | Influence of urban land development and subsequent soil rehabilitation on soil aggregates, carbon, and hydraulic conductivity[END_REF]. In our case, the effect may disappear when Anthrosols and Technosols are rich in organic matter (SOC > 3.4 % and N > 2.4 g kg -1 ) and their current land-use is older than 13 years. SOC and N values above a threshold of 3.4 % and 2.4 g kg -1 respectively may not induce further positive effect on Kfs, but it is observed that amendments increase Kfs in most of cases [START_REF] Kranz | The effects of compost incorporation on soil physical properties in urban soils -A concise review[END_REF]. years, but it was twice as high as the porosity after 40 years (70 % v/v). Such an increase in porosity increased the initial Kfs from 11.7 to 66.7 µm s -1 . Thus, the evolution of porosity may explain why Anthrosols and Technosols with an older land-use history have larger Kfs. Lastly, the way coarse materials are incorporated into the soil or deposited at the surface, the lack of ploughing or amendments may reduce Kfs [START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF]. Therefore, the assessment of the infiltration on vegetated urban soils requires a detailed study of the historical background.

The type of vegetation cover is not relevant for run-off estimation

The type of vegetation cover might modify stormwater infiltration in Anthrosols and Technosols. We hypothesized that soils presenting a mixture of trees, shrubs and herbaceous plants would have high Kfs values, as found in [START_REF] Wang | Effects of Urbanization, Soil Property and Vegetation Configuration on Soil Infiltration of Urban Forest in Changchun, Northeast China[END_REF]. Although it was easy to establish categories of vegetation cover, it proved to be difficult to identify them in the field because the vegetation cover was heterogeneous and often had several configurations and densities for the same plot.

As a result, the type of vegetation was not a determining factor for the stormwater infiltration, as it was also found by [START_REF] Yang | Water infiltration in urban soils and its effects on the quantity and quality of runoff[END_REF]. In addition, the type of vegetation cover does not have a significant role in run-off regulation in natural Mediterranean soils, because bare soils are the main generators of run-off [START_REF] Mayor | Factors and interactions controlling infiltration, runoff, and soil loss at the microscale in a patchy Mediterranean semiarid landscape[END_REF].

In order to better understand how to mitigate the risk of run-off flood in urban areas, it would be necessary to study Kfs on bare Anthrosols and Technosols to check the data against natural soils. Bare urban soils may likely produce significant run-off in areas subjected to episodes of heavy precipitation. To prevent run-off, the vegetation cover should be maintained in these areas. However, spatial planning should adapt the management to the properties of Technosols because the low available water content can expose plants to hydric stress [START_REF] Huot | Lysimeter monitoring as assessment of the potential for revegetation to manage former iron industry settling ponds[END_REF][START_REF] Urban | Two Different Approaches to Improve Growing Conditions for Trees Comparing Silva Cells and Structural Soil[END_REF].

Conclusion

Our results showed that vegetated urban soils can infiltrate enough stormwater to mitigate run-off in urban areas with episodes of heavy precipitation, such as Marseilles (France). The estimation of potential levels of run-off production in Anthrosols and Technosols based on the soil depth and the field-saturated hydraulic conductivity (Kfs) enabled soils to be discriminated according to their land-use history. Organic matter appeared to have a positive effect on Kfs in the early stage of T soils formation and hence on reducing their potential run-off. Urban soils are a key asset to mitigate the risk of run-off flood in urban areas, but are usually neglected in planning documents, while we show that they may have a low level of run-off production via a large Kfs (> 10 µm s -1 ) and a significant depth (> 0.5 m).

Therefore, a sound planning approach to mitigate run-off should consider Anthrosols and Technosols, along with other strategies (e.g. stormwater basins).

This implies characterizing these soils for their physico-chemical and physical properties, including their infiltration, but also assessing their land-use history. In doing so, it would contribute to the new European strategy to minimize run-off in areas prone to extreme episodes of precipitation (European Commission, 2015).

Although quantifying the run-off was beyond the scope of this paper, the minimization of soil sealing appears to be a sustainable strategy to mitigate run-off.

Further research is needed, however, to refine the estimate of field-saturated soil hydraulic conductivity (Kfs) in Anthrosols and Technosols of urban areas, since the method used to measure Kfs can affect the estimation of the run-off.

  Local d'Urbanisme (PLU), Schéma de cohérence territoriale (SCOT) -and flood risk prevention strategies -Stratégie Locale de Gestion de Risques d'Inondation (SLGRI), Plan de Prévention des Risques d'Inondation (PPRI) -do not yet consider that the hydraulic properties of urban soils may influence stormwater infiltration.

  photointerpretation of the 2011 aerial images, the most recent images at the

4. 4 .

 4 Impact of land-use history on infiltration in urban soilsThe land-use history affects the Kfs of Anthrosols and Technosols, and thus the run-off production. According to our results, Kfs was low in recent urban soils with a current vegetated land-use < 13-year-old. The influence of the historical background could be due to several factors, such as the intensity of earthworks, the structure evolution of Technosols, the way in which construction debris was deposited and soil was managed. In our set of soils, the low Kfs corresponded to Anthrosols and Technosols that had experienced intense earthworks during dismantling of the tile factory and quarries. Similarly,[START_REF] Shuster | Residential demolition and its impact on vacant lot hydrology: Implications for the management of stormwater and sewer system overflows[END_REF] found, in Technosols subjected to demolition, that < 14-15-year-old soils had a lower infiltration than > 14-15-year-old soils, since recent Technosols had been affected by more intense earthworks.The structure of vegetated urban soils may evolve with the age of the current landuse. The structure of one-year-old constructed Technosols and Anthrosols can be dominated by gravity and compaction and, is more influenced by roots and soil organisms from 2-4 years on[START_REF] Jangorzo | Evolution of the pore structure of constructed Technosols during early pedogenesis quantified by image analysis[END_REF][START_REF] Bouzouidja | Green roof ageing or Isolatic Technosol's pedogenesis?[END_REF]. Seré et al. (2012) tested a constructed Technosol in lysimeters for 3 years. They observed an early process of pedogenesis, which created sub-horizons and improved the structure by the formation of micro and macro-aggregates.[START_REF] Scalenghe | The First Forty Years of a Technosol[END_REF] studied three Technosols on a slope embankment of 0, 4 and 40 years. They also observed the formation of horizons and described in detail the evolution of porosity in a Technosol. The initial total porosity was reduced after 4

Table 3

 3 Urban soils characteristics. Descriptive statistics for the first horizon of Soil Organic Carbon content (SOC), Total Kjeldahl nitrogen (N), and Olsen phosphorus (P); Descriptive statistics of the infiltration rate (Infiltration), field-saturated hydraulic conductivity (Kfs), and total depth of the soil profile (Depth); soils with estimated low run-off production (Low run-off); for (T) plots that lay on former clay quarries, tile factories; for (NT) plots that did not lie on former clay quarries or tile factories (NT), and (C) shrubland used as control soil.

		SOC	N	P	Infiltration	Kfs	Depth	Low run-off
		(%)	(g kg -1 )	(g kg -1 )	(mm h -1 )	(µm s -1 )	(m)	(%)
	T (n=14)							84.6
	Geometric Mean ± SD	2.7 ± 2.8	2.0 ± 1.8 0.05 ± 0.08 117.6 ± 546.7 12.3 ± 51.0 0.59 ± 0.28
	Median	3	2.2	0.06	54.1	6.0	0.68
	Min	0.2	0.3	0.01	14.4	1.5	0.20
	Max	10.8	7.3	0.29	1705.3	154.7	1.00
	NT (n=13)							14.3
	Geometric Mean ± SD	5.4 ± 1.6	3.0 ± 1.4 0.33 ± 0.19 301.3 ± 190.2 26.7 ± 16.0 0.71 ± 0.23
	Median	5.3	2.8	0.29	333.3	29.5	0.70
	Min	3.4	2.0	0.14	71.7	7.3	0.45
	Max	9.2	7.0	0.80	707.1	57.5	1.10
	C (n=1)
	Value	3.9	2.5	0.01	2334.9	150.1	0.40

Table 4 :

 4 Pearson Correlation values between Field-saturated soil hydraulic conductivity (Kfs) and Soil Organic Carbon (SOC), Total Kjeldahl N (N), Olsen P (P), clay content (Clay), silt content (Silt), sand content (Sand), and Penetration resistance in the first 0.15 m (Qd). The Pearson correlation (r) values are shown for: (1) plots that lay on former clay quarries, tile factories (T), and(2) plots that did not lie on former clay quarries or tile factories (NT). Pearson coefficients statistically significant (P-values < 0.050) are noted as *. Observed trends between variables which were not statistically significant (0.051 < Pvalue < 0.069) are noted as 1 .

	Kfs	SOC	N	P	Clay	Silt	Sand	Qd
	T n=14	0.50 1	0.55*	-0.16	0.50 1	0.15	-0.37	-0.41
	NT n=13	0.46	0.54 1	0.12	-0.21	-0.15	0.18	-0.35

de Riaux-Estaque and the inhabitants of the 16 th district of Marseilles. We thank Jeanne Maréchal, Adèle Gomez and Bernard Angeletti for their help in the field and lab work. We thank Daniel Vázquez-Tarrio and Jeanne Maréchal for their comments on the permeameter design. We also thank an anonymous reviewer for his/her very helpful comments. Table 1 Suggested

Tables


B Leistedt

A Rassat

A Réfrégier

J.-L Starck

3DEX: a code for fast spherical Fourier-Bessel decomposition of 3D surveys

Keywords: Cosmology, HEALPix, spherical Fourier-Bessel decomposition, spherical signal processing

3DEX: a code for fast spherical Fourier-Bessel decomposition of 3D surveys

Introduction

In the last few decades, cosmology has become a datadriven field, where high-precision measurements of the cosmic microwave background (CMB, e.g., Larson et al. 2011), weak lensing (e.g., Schrabback et al. 2010) and galaxy surveys (e.g., Percival et al. 2007 b) have permitted the establishment of a standard cosmological model in which the Universe is composed of 4% baryons, 22% dark matter and 74% dark energy. Some major questions remain, the nature of dark matter and dark energy in particular is still not understood. Similarly, the initial conditions of the Universe are yet to be established and alternative models of gravity are still to be tested in comparison with Einstein's general relativity.

New surveys are underway with these science objectives, e.g. Planck for the CMB (The Planck Collaboration 2006), DES (Dark Energy Survey, Annis et al. 2005), BOSS (Baryon Oscillation Spectroscopic Survey, Schlegel et al. 2007), LSST (Large Synoptic Survey Telescope, Tyson & LSST 2004) and Euclid (Laureijs et al. 2011;Refregier et al. 2010) for weak lensing and the study of large-scale structure with galaxy surveys. In order to be beneficial, cosmological boris.leistedt@gmail.com anais.rassat@epfl.ch studies of these surveys need to use high-precision statistical methods, such as a full 3D analysis on the sky where all-sky 3D surveys are available.

Several tools have been developed to analyse data on the sphere, which is required for a 2D spherical harmonic CMB analysis (Crittenden & Turok 1998;Crittenden 2000;Górski et al. 2002Górski et al. , 2005;;Doroshkevich et al. 2008). Weak lensing and galaxy survey data can also be analysed tomographically (i.e. in 2D slices), but unlike for the CMB, a full 3D spherical Fourier-Bessel analysis can also be sought (Fisher et al. 1995;Heavens & Taylor 1995;Heavens 2003;Castro et al. 2005;Rassat & Refregier 2011). Previous 3D data analyses were on relatively small data sets (Fisher et al. 1995;Heavens & Taylor 1995;Erdogdu (b) et al. 2006;Erdogdu (a) et al. 2006), but future surveys like Euclid and LSST will provide surveys with billions of galaxies, making previous methods for calculating the 3D spectra unfeasibly time-consuming.

In Section 2.1, we present the theory behind the 3D Fourier-Bessel decomposition for infinite and finite continuous fields as well as the usual method for a discrete survey (e.g. galaxy survey). In section 2.2, we present two additional equivalent formulations of the spherical Fourier-Bessel decomposition, one of which is central to the 3DEX code. In Section 3, we compare the accuracy and calcula-tion time for the usual method used for calculating Fourier-Bessel coefficients and methods with the 3DEX code presented in this paper. In Section 4, we describe the 3DEX library and give examples of how to use it. In Section 5 we present our conclusions. We also include an appendix, where we discuss the subtleties of the Fourier-Bessel normalisation.

Theory

The spherical Fourier-Bessel decomposition

In observational cosmology, spherical coordinates (where the observer is at the origin) are a natural choice for the analysis of cosmological fields. In this system of coordinates, eigenfunctions of the Laplacian operator are products of spherical Bessel functions and spherical harmonics, i.e. functions j (kr)Y m with eigenvalues -k 2 . For an homogeneous three-dimensional field f (r) = f (r, θ, φ) in a flat geometry, the spherical Fourier-Bessel decomposition (Fisher et al. 1995;Heavens 2003;Castro et al. 2005) is

f (r, θ, φ) = 2 π dk m f m (k)kj (kr)Y m (θ, φ), (1)
with the inverse relation

f m (k) = 2 π d 3 r f (r, θ, φ)kj (kr)Y * m (θ, φ). (2)
Note that this decomposition uses the same notation as Rassat &Refregier (2011) andCastro et al. (2005), which is slightly different from the one used in Lanusse, Rassat, & Starck (2011). The coefficients may be used to calculate the 3D power spectrum C( , k), defined by

f m (k)f * m (k ) = C(l, k)δ D (k -k )δ δ mm , (3)
a naïve estimator of which is

C (k) = 1 2l + 1 m |f m (k)| 2 . ( 4
)
This can be seen as an extension of the usual 2D power spectrum f m f * m = C l δ δ mm . The latter arises from the spherical harmonic transform of a 2D field given on the sphere f (θ, φ) = m f m Y m (θ, φ). In practice, surveys will only cover a finite amount of volume, limiting the analysis to a sphere of radius R. These boundary conditions lead to a discrete spectrum {k n }, which is detailed in the appendices. In this paper, we assumed as a boundary condition that f vanishes at r = R. The spherical Fourier-Bessel decomposition becomes (Erdogdu (b) et al. 2006;Fisher et al. 1995)

f (r, θ, φ) = mn κ n f m (k n )k n j (k n r)Y m (θ, φ), (5)
which is exact if the ranges of ,m and n are infinite. The Fourier-Bessel coefficients are denoted by f mn = f m (k n ), and κ n is the normalisation constant (see appendices for more details).

In various applications, though, the continuous field f cannot be directly observed. This is notably the case in cosmology where galaxy surveys give indirect information about the underlying matter density field through their spacial positions. Note that these tracers are subject to various distortions and non-linearities, but these are not the purpose of this work. In this work we only consider linear or quasi-linear scales ( < 50, k < 0.2hMpc -1 ).

If the only information about the field f is a list of coordinates r p = (r p , θ p , φ p ) with p = 1, . . . , N (where N is the number of galaxies in the latter example), the survey may be considered as a superposition of 3D Dirac deltas and each coefficient f mn can simply be estimated with a sum (Heavens & Taylor 1995;Fisher et al. 1995;Erdogdu (b) et al. 2006;Abramo et al. 2010)

f (r) = N p=1 δ (3) (r -r p ), (6)
f mn = N p=1 k n j (k n r p )Y * m (θ p , φ p ). ( 7
)
This formulation has been used for the analysis of shallow galaxy surveys such as the IRAS 1.2mJ survey (∼ 6k galaxies, Strauss et al. 1992;Fisher et al. 1995;Heavens & Taylor 1995), and the 2MRS survey (2MASS Redshift Survey, ∼ 45k galaxies, Huchra et al. 2011;Erdogdu (b) et al. 2006;Erdogdu (a) et al. 2006). Since the time to calculate equation 7 is proportional to N n max ( max + 1) 2 /2, Equation 7 will become highly time-consuming when applied to larger surveys or when precise decomposition is required (large n max and max ).

Three equivalent formulations

In spherical coordinates, since 3D space can be viewed as an infinite series of closed shells Ω(r), the spherical Fourier-Bessel decomposition may also arise from repeated 2D spherical harmonic transforms to which spherical Bessel transforms are applied (Abramo et al. 2010). Formally, the field f given on each shell Ω(r) is first expanded into spherical harmonics

f (r, θ, φ) = m f m (r)Y m (θ, φ), (8)
for which the inversion formula gives harmonics coefficients f m (r) depending on the radius r

f m (r) = Ω(r) dΩ f (r, θ, φ)Y * m (θ, φ). (9)
It is then possible to perform a spherical Bessel transform

f m (r) = 2 π dk f m (k)kj (kr), (10)
leading to the final Fourier-Bessel coefficients f m (k)

f m (k) = 2 π dr r 2 f m (r)kj (kr). ( 11
)
This formulation hence extends the notion of 2D spherical harmonics to three-dimensional fields.

It is also possible to conceive the reverse approach, i.e. to perform the spherical Bessel transform first and subsequently expand the resulting coefficients into spherical harmonics. Formally, the -th order spherical Bessel transform of f (similar to its Hankel transform) is

f (r, θ, φ) = 2 π dk f (k, θ, φ)kj (kr), (12)
for which the inversion formula gives

f (k, θ, φ) = 2 π dr r 2 f (r, θ, φ)kj (kr). ( 13
)
The result is then expanded into spherical harmonics but with an unusual formulation since f (k, θ, φ) and Y m (θ, φ) (as well as the basis functions j (kr) and Y m (θ, φ)) have now the parameter in common:

f (k, θ, φ) = m f m (k)Y m (θ, φ). ( 14
)
Again, using the inversion formula, we obtain the Fourier-Bessel coefficients

f m (k) = Ω dΩ f (k, θ, φ)Y * m (θ, φ). ( 15
)
Due to the closed domains of shells Ω(r) and thus the relative independence of angular and radial dimensions, the raw (equations 1 and 2), the forward (denoted by SHB for spherical-Harmonic-Bessel, equations 8 to 11) and the reverse (denoted by SBH for spherical-Bessel-Harmonic, equations 12 to 15) methods are equivalent formulations of the spherical Fourier-Bessel decomposition of any threedimensional field f (r, θ, φ). This is summarised in the following schematic description of each method:

RAW : f (r) three-dimensional integration ------------------→ f m (k) SHB : f (r) SHT ---→ f m (r) SBT ---→ f m (k) (16) SBH : f (r) SBT ---→ f (k, θ, φ) SHT ---→ f m (k).
Note that this section is related to the ideal case R = ∞, but all equations can be straightforwardly rewritten for a finite R by replacing k by k ln , bounding each integral and adapting normalisation. The formulas arising from this adaptation are used in the next sections.

Estimating Fourier-Bessel coefficients from a real survey

Although the three approaches described in 2.2 are theoretically equivalent, their estimates and numerical implementations take different forms.

Forward approach (SHB)

Estimating the f mn coefficients using the forward method naturally requires the radial dimension to be discretised. Indeed, the first step is to compute the spherical harmonic transform on a set of shells located at radial values r 1 , . . . , r N layers . In each layer, the coefficients f m (r i ) are estimated. Although it is possible to perform a raw estimate for the later harmonics transform, it is often advisable to use a robust 2D discretisation scheme (of N pix (i) pixels for the i-th shell) and to take advantage of the related high-performance algorithms. Angular space is hence discretised into nodes (r i , θ p , φ p ) = (r i , γ q ) and the field is approximated on each node, giving f (r i , γ p ). The spherical harmonic decomposition in the i-th shell becomes

f m (r i ) = N pix(i) p=1 f (r i , γ p )Y * m (γ p ), (17)
and the final coefficients are obtained by performing the following spherical Bessel decomposition:

f mn = N layers i=1 f m (r i )k n j (k n r i ). ( 18
)
With this method, radial and angular spaces are discretised and both transforms are approximated.

Reverse approach (SBH)

For the reverse approach, a 2D scheme on the sphere was required as well. As previously, this scheme defines a set of N pix zones (pixels) related to angular nodes γ q . If G q denotes the points of the survey located in the solid angle corresponding to the q-th zone of the scheme, we perform the spherical Bessel Transform (raw estimate) in each zone

f n (γ q ) = f (k n , γ q ) = p∈Gq k n j (k n r p ), ( 19
)
and each of these intermediate maps is decomposed into spherical harmonic (spherical Harmonics Transform) which gives the Fourier-Bessel coefficients

f mn = N pix q=1 f n (γ q )Y * m (γ q ). ( 20
)
With the reverse method, one can avoid to discretise radial space. Moreover, this one-shell pixelisation of the sky (thus based on physical solid angles) allows for a natural treatment of radial distortions (redshift, relativistic) and masking effects. Using multiple resolutions at different radial values, as would be possible with the forward method, is much more questionable. The SHB method also proves to be a powerful tool for weighting the data prior to estimating the power spectrum. For instance, in Tadros et al. (1999) used a fiducial power spectrum to derive an optimal weighting operation. This operation is quite complex when using the raw Fourier-Bessel approach, whereas the SHB formulation naturally handles the dependence on k of the weighting function.

The three methods to estimate the spherical Fourier-Bessel decomposition can therefore also be expressed for a discrete 3D survey, summarised schematically below:

RAW : {r p } Raw sum, best estimate of FB coefficients --------------------------→ f mn SHB : {r p } Approx SHT --------→ f m (r i ) Approx SBT --------→ f mn SBH : {r p } Exact SBT -------→ f n (γ p ) Approx SHT --------→ f mn .
Note that in practice, the range of (l, m, n) is finite, which introduces an additional approximation. Here, and n are restricted to [0, max ] and [1, n max ] respectively. Given , m goes from -to .

Method comparison

Complexity, accuracy and discretisation grids

For a survey that probes a field by N discrete points, the raw method is the natural estimate of the Fourier-Bessel coefficients. However, since each point contributes to the calculation of every coefficient f mn (∀ l, m, n), computation time is proportional to N • n max ( max + 1) 2 /2, which can be highly problematic for large surveys.

In the forward method, the repeated spherical harmonic transforms take advantage of tesselation schemes and high-performance algorithms such as those provided by HEALPix [Górski et al. (2005)], IGLOO [Crittenden & Turok (1998)] or GLESP [Doroshkevich et al. (2008)]. Roughly speaking, the number of nodes to be considered is reduced from N to N pix , and the use of fast spherical harmonic transforms on these schemes significantly decreases computation time.

However, this approach requires the three-dimensional space to be divided into shells Ω(r i ). Both radial and angular dimensions are discretised, and the survey is approximated on an actual 3D grid. In practice, this approximation deteriorates the accuracy of the estimated Fourier-Bessel coefficients. Furthermore, designing a meaningful radial discretisation is a difficult task. For equal-area pixelisations, the area of each pixel on the i-th shell is 4πr 2 i /N pix (i). With HEALPix, the n side angular parameter may only be increased by a factor 2, which changes the number of pixels by a factor of 4 (since N pix (i) = 12n side (i) 2 ). This means that pixel areas cannot be stabilised for subsequent shells as r increases. Consequently, it is not possible to adapt a resolution to obtain a 3D scheme with equal-volume voxels. Extending 2D schemes to 3D is difficult and may even require a novel formalism for an equal-voxel 3D grid.

In the reverse approach, though, the use of angular 2D schemes is possible, but radial space does not need to be discretised. The points of the survey are grouped according to angular zones instead of being approximated on a 3D grid. An estimate of the spherical Bessel transform is computed in every solid angle, and the result is then expanded in spherical harmonics on the 2D spherical grid. In the final account, this method naturally leads to more accurate coefficients than the forward method and also takes advantage of high-performance 2D schemes. For these reasons, we focus on the reverse approach and its implementation, using HEALPix for the angular transform.

Finally, for both forward and reverse methods the spherical harmonics discretised basis (coefficients Y m (γ p )) may be fully pre-computed and stored in external files. This is a particularly useful feature (incompatible with the original formulation of spherical Fourier-Bessel), which significantly eases and speeds up the use of these methods.

Speed and accuracy of the Reverse Method

To test the accuracy and speed of the reverse method compared to the raw method, we considered the high-resolution full-sky Horizon simulation (Teyssier et al. 2009). Horizon is a N-body simulation covering a 2h -1 Gpc periodic box using 70 billion dark matter particles using a WMAP3 cosmology (Spergel et al. 2007). A derived halo catalogue is available, which we used to calculate f mn and C (k n ) values using both methods (raw and reverse). Since we are interested only in comparing the speed of each method, we simply considered each halo to have equal weight.

We performed the raw and the reverse estimates on three 'surveys' of N = 4.2 × 10 5 , 3.1 × 10 6 and 1 × 10 7 halos, which correspond to three different depths (z max = 0.1, 0.2 and 0.3 respectively) in the Horizon simulation. The HEALPix angular parameter is given by n side .

The results of the accuracy and speed tests are given in Table 1. The third (fourth) column gives the percentage f coefficients for which the relative accuracy (f mn ) ( (C n )) is lower than 0.3% for given values of n side and N . We considered the intervals (l, n) ∈ ([0, 20], [1,20]) and (l, n) ∈ ([21, 50], [21,50]) separately, since the estimation of higher coefficients depends more on the value of n side . We also compared computation times of the two methods by observing the ratio T = t reverse /t raw . Given a survey and a method, computation time denotes the CPU time required to compute the k ln 's (from the Bessel functions) and the final coefficients f mn without using pre-computed quantities. Note that we performed this analysis by distributing the calculations on five machines and simply adding the individual contributions to computation time since our method is linear with survey size. With the reverse method, though, the roots of the Bessel functions as well as the spherical harmonics may be pre-computed and stored in external files, which decreases computation time and complexity when working with 3DEX. The reverse method is about an order of magnitude faster than the raw method, but this depends on the choice of n side . For n side = 1024 almost all f mn coefficients in the range [0, 20] (for and n) have relative error below 1%, and 90% have it below 0.3%, whereas over 99% of C( , k n ) coefficients are accurate to < 0.3%. In the range [20,50], the accuracy is somewhat degraded due to the extension of the HEALPix formalism to 3D surveys. Indeed, for data distributed on the sphere, 3D space is very sparse even for large surveys. Increasing n side to 1024 or 2048 strongly improves the accuracy for higher orders . Note that comparisons for > 50 are limited by the amount of time the raw method takes. Fig. 1. Speed results (raw and reverse methods) for increasing summation limits max = n max , for a survey of N = 9.7×10 6 halos. Dashed lines are the fitted complexity curves. The reverse formulation is suitable for pre-calculations and external storage of the spherical harmonics, which was not performed here but enables for additional speed improvements.

N n side r (f mn ) < 0.3% r (C n ) < 0.3% T [0,20] / [
Figures 1 and2 show the time taken for calculations as a function of max and n max (Figure 1), and as a function of number of halos (Figure 2). The boxes correspond to the raw method, the circles and diamonds to the reverse method with n side = 512, 1024 respectively. The dashed line corresponds to the general rule that the raw method scales as N n max ( max + 1)

2 /2, whereas the points are all estimated from calculations. With max = n max = 100 and N = 9.7e6, the raw decomposition took a few days of calculations, whereas the reverse method only took 12 hours. In our formalism, k max is determined by the choice of R for the boundary condition and by the band-limit n max for spherical Bessel coefficients. For each multipole we have k max = k nmax = q nmax /R where q nmax is the n maxthe root of the -th spherical Bessel function. Because R is usually imposed by the problem or the data, one must increase n max to probe smaller radial scales. In fact, a reasonable approximation (or even a simple plot) shows that q n ≈ ( + 3n). This observation enabled to predict which radial scales are probed and how computation time scales with k max , given that we provide the complexity for max and n max .

One of the main advantages of the reverse method is that it is naturally suited to parallel computing it uses HEALPix fast spherical Harmonics Transform routines. All previous tests were performed on a recent computer (i7 processor, 8Go RAM) and take advantage of OpenMP (with four threads). More advanced computing means (larger RAM and more processors) significantly decrease calculation time. For example, max = n max = 128 with n side = 2048 took about an hour with 128 cores and 512Go RAM, whereas computation time for the raw method was estimated to several days on the same machine. Note that the raw method is also suited to parallelisation: galaxies may be treated separately by different threads. In all experiments, we took advantage of this property and performed both raw and reverse decompositions with four threads to perform relevant comparisons between the two.

In terms of the power spectrum, figures 3 and 4 show the relative error between the raw and the reverse methods both in mode-mode space ( -n) and in mode-scale space ( -k n ). For this comparison we decomposed a survey of N = 4.2×10 5 halos with z max = 0.1. Figures in mode- In view of the -k n space, we see that no fluctuations are observed along the k axis up to k = 0.03hMpc -1 . In this range, fluctuations occur in space, which are accurately probed with n side = 512 until = 25 but naturally require a more precise scheme for > 25, k > 0.03hMpc -1 (smaller scales in physical space). In conclusion, parameter n side (as well as R) must be chosen depending on the scales one wishes to probe. Figures 3 and4 provide accuracy results that are complementary to Table 1.

The 3DEX library

The 3DEX library requires the HEALPix package (v2.12 or later) and the CFITSIO library. 3DEX may either by installed with an HEALPix-like procedure (configure and make commands) or using CMake. The Fortran modules, the 3DEX dynamic library and the related executables will be created in the relevant directories (see README file for more information).

In addition to the numerical procedures required to compute Fourier-Bessel coefficients, various other routines are provided in the library, such as those converting redshift to comoving distance, computing spherical Bessel functions and their zeros, reading and writing 3D structures (f lmn and C ln ), or reconstructing radial maps from Fourier-Bessel coefficients.

Three executable programmes are generated:

-survey2almn performs the spherical Fourier-Bessel decomposition (reverse method) of a discrete survey with input parameters l max , n max , r and n side . Outputs are the f lmn coefficients and the power spectrum. -survey2almn interactive is very similar to the previous programme, but converts redshift values into comoving distance before performing the spherical Fourier-Bessel decomposition. In particular, the routine takes a .txt file as input, taking into account parameters on the cosmology and on the decomposition. -almn2rmap extracts the f lmn coefficients from a FITS file and reconstructs the field (HEALPix map) at a given radius. Inputs are the resolution, the radius and summation limits l max and n max , which allows one to reconstruct several maps at different scales and resolutions.

The corresponding calls are given by the examples below.

> survey2almn survey thetaphir.dat almn.fits cln.fits 20 20 256 2000.0, where survey thetaphir.dat is a survey with columns representing θ, φ, r, and the keywords correspond to values of , n, n side , R. The output is both the coefficient values (almn.fis) and the Fourier-Bessel spectrum (cln.fits).

> survey2almn interactive parameters.txt, where parameters.txt is an external file containing input parameters for the survey and the cosmology (which allows for more flexible use). Finally, for the map reconstruction, we can use: > almn2rmap almn.fits map.fits 400.0 256 10 10 2000.0, where the keywords correspond to r max , n side , , n, R.

Conclusion

High-precision cosmology from galaxy and weak lensing surveys will require the analysis of 3D data in spherical coordinates, a situation for which spherical Fourier-Bessel decomposition is most suited. Current methods will be inadequate for future planned cosmological surveys, which will provide for example galaxy surveys with billions of galaxies, compared to millions today.

We have reviewed the forward or SHB (spherical Harmonic-Bessel) formalism of the spherical Fourier-Bessel decomposition which first calculates the tangential, then the radial decomposition. We also introduced the reverse or SBH (spherical Bessel Harmonic) formalism that inverses this order. Only the latter approach can take advantage of existing fast codes for the calculation of tangential modes. (To do the same, the former would require a a new voxelisation scheme.)

We presented a public code 3DEX (3D EXpansions) for the fast calculation of Fourier-Bessel coefficients and spectra, which uses the HEALPix pixelisation scheme for calculating the tangential modes. The 3DEX code is based on the reverse/SBH formulation of the Fourier-Bessel decomposition.

We tested the 3DEX code on linear and quasi-linear scales ( < 50 and k n < 0.2hMpc -1 ) using the Horizon halo simulation for redshifts z < 0.3. For n side = 1024 the 3DEX method for calculating the power spectrum C( , k) is accurate to 0.3% on these scales.

For surveys with < 10 million galaxies, computation time is reduced by a factor 4-12 depending on the desired scales and accuracy. For larger surveys the gain in time will be even greater. Finally, the use of the 3DEX code is not restricted to cosmological calculations, and can be used in any other discipline that requires a spherical Fourier-Bessel analysis of 3D data.

Fig. 2 .

 2 Fig. 2. Speed results (raw and reverse methods) for increasing survey size, for max = n max = 30. Dashed lines are the fitted complexity curves.

  mode space naturally differ according to the choice of the boundary R because the latter determines the discrete radial scale spectrum {k n }, and hence mode n computed with two different R's corresponds to different k-scales. When comparing the results from R = 1000 and R = 2000 in mode-scale space, we observe that the boundary condition fixes the explored scales. The left column is thus complementary to the right column to explore higher values of k. Although figure 3 gives information about the final coefficients, figure 4 is hence more appropriate to see which scales are probed and with what accuracy.

Fig. 3 .

 3 Fig. 3. Relative error on the power spectrum in mode-mode space C(l, n). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing n side = 512, 1024 (rows) and R = 1000, 2000 (columns). Only a few zones (white spots) are outside the range [-0.3%, +0.3%].

Fig. 4 .

 4 Fig. 4. Relative error on the power spectrum in mode-scale space C(l, k) (k is in hMpc -1 ). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing n side = 512, 1024 (rows) and R = 1000, 2000 (columns). Only a few zones (white spots) are outside the range [-0.3%, +0.3%].

Table 1 .

 1 Estimation of Fourier-Bessel coefficients: comparison of the new method, the reverse formulation (equations 19 and 20 using HEALPix discretisation) with the original, raw formulation (equation 7). The third (fourth) column gives the percentage f coefficients for which the relative accuracy (f mn ) ( (C n )) is lower than 0.3% for given values of n side and N . T is the ratio of elapsed times of the two methods.

	21,50]	[0,20] / [21,50]

Acknowledgements. The authors are grateful to Ofer Lahav, Pirin Erdogdu and François Lanusse for useful discussions regarding the spherical Fourier-Bessel theory, as well as to Romain Teyssier and Nicolas Clerc for access and help using the Horizon simulation. The authors thank Yves Revaz for computational help at EPFL. The 3DEX library uses Healpix software (Górski et al. 2002(Górski et al. , 2005)). This work is supported by the European Research Council grant SparseAstro (ERC-228261). This research is in part supported by the Swiss National Science Foundation (SNSF).

Appendix A: Normalisation and discrete radial spectrum

The basis functions kj l (kr)Y lm (θ, φ) form a set of eigenfunctions of the Laplacian operator in spherical coordinates. In particular, these functions are orthonormalised in the continuous case, thanks to the orthogonality relation (Baddour 2010) where δ K is Kronecker's delta notation and δ D Dirac's function.

A common approach to simplify the problem is to assume some boundary conditions for the field f . Different conditions have been explored in the literature (Fisher et al. 1995;Heavens & Taylor 1995), including potential or gradient continuity. In this paper, we used a condition that derives from the classical formulation of the discrete spherical Bessel transform: space is assumed to be finite and limited to a sphere of radius R. In this case, the spherical Bessel functions are not normalised and the boundary effect leads to a discrete spectrum {k ln }. The Fourier-Bessel coefficients become a set f lmn = f lm (k ln ), and the complete description of the field in the so-called Fourier-Bessel basis (Binney & Quinn 1991) is summarised in equation 5.

As a consequence, a natural choice for the boundary condition is to impose the field to vanish at r = R (Abramo et al. 2010), which constrains the Bessel functions and generates the radial spectrum {k ln } such that, for all l and n,

If q ln denotes the n-th root of j l (z), the closure relation of the Bessel basis is

The discrete spectrum is thus fixed by the zeros of the spherical Bessel functions. We obtain the normalisation coefficients κ ln (Fisher et al. 1995)

which are used for field reconstruction (equation 5).

Other approaches are possible to tackle boundary conditions in radial space, notably those imposing potential continuity at r = R (Fisher et al. 1995). Then, the discrete spectrum k ln is such that

and normalisation constraint becomes

Appendix B: Angular masks 3DEX takes into account optional angular masks under the form of either an equatorial cut or an input all-sky FITS map.

In the first case, supplying θ cut defines the latitude (in degrees) of a straight symmetric cut around the equator. Pixels located within that cut (l = cos(θ cut )) are ignored.

In the second case, the supplied mask must be an HEALPix map (ring ordering) of N pix pixels at resolution n side (which must be identical to Fourier-Bessel resolution parameter) w(γ q ) q=1,...,Npix .

(B.1)

In the forward method, the first step is to apply the mask to each discrete shell before performing the spherical Harmonics Hence for the ith shell, field f is weighted by the mask at each pixel

The Fourier-Bessel coefficients are obtained after performing SH and SB transforms.

In the reverse method, the first step is still the spherical Bessel Transform, which gives a set of n max HEALPix maps f (k ln , γ q ). The mask is then applied to each of these maps f (k ln , γ q ) = w(γ q )f (k ln , γ q ). (B.3) and the modified spherical Harmonics Transform gives the final f lmn coefficients.