P.-A Reninger

G Martelet

J Perrin

J Deparis

Y Chen

Slopes of an airborne electromagnetic resistivity model interpolated jointly with borehole data for 3D geological modelling

Keywords: Modelling, Airborne electromagnetic, Resistivity sounding

We investigate a novel way to introduce resistivity models deriving from airborne electromagnetic surveys into regional geological modelling. Standard geometrical geological modelling can be strengthened using geophysical data. Here, we propose to extract information contained in a resistivity model in the form of local slopes that constrain the modelling of geological interfaces. The proposed method is illustrated on an airborne electromagnetic survey conducted in the region of Courtenay in France. First, a resistivity contrast corresponding to the clay/chalk interface was interpreted confronting the electromagnetic soundings to boreholes. Slopes were then sampled on this geophysical model and jointly interpolated with the clay/chalk interface documented in boreholes using an implicit 3D potential-field method. In order to evaluate this new joint geophysical-geological model, its accuracy was compared with that of both pure geological and pure geophysical models for various borehole configurations. The proposed joint modelling yields the most accurate clay/chalk interface whatever the number and location of boreholes taken into account for modelling and validation. Compared with standard geological modelling, the approach introduces in between boreholes geometrical information derived from geophysical results. Compared with conventional resistivity interpretation of the geophysical model, it reduces drift effects and honours the boreholes. The method therefore improves what is commonly obtained with geological or geophysical data separately, making it very attractive for robust 3D geological modelling of the subsurface.

I N T R O D U C T I O N

Three-dimensional (3D) geological modelling, often used in the oil industry, is increasingly being adapted in various domains [START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method[END_REF]. It is a useful tool to quantitatively represent the extension of geological structures [START_REF] Turner | Challenges and trends for geological modelling and visualization[END_REF]. It commonly results from the interpolation of the information available in boreholes, cross-sections, and/or geological maps [START_REF] Kaufmann | 3D geological modelling from boreholes, cross-sections and geological maps, application over former natural gas storages in coal mines[END_REF]. Recent geostatistical developments have benefitted from geological * E-mail: pa.reninger@brgm.fr modelling since orientation and dip of geological formations observed in the field can now be considered to constrain 3D geometries using the potential-field method [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation[END_REF][START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules[END_REF]. However, often, geological data are scarce and scattered, inducing large uncertainties on 3D geological models. Geophysics is a good candidate to complement geological data at depth since it provides supplementary information to surficial data. In previous studies, with the use of seismic data or the inversion of gravity, magnetic data have been proved to be useful to reduce uncertainty at depth and thus to define more robust geological models (e.g., [START_REF] Martelet | Integrated 3D geophysical and geological modelling of the Hercynian Suture Zone in the Champtoceaux area (South Brittany, France)[END_REF][START_REF] Guillen | Geological modelling from field data and geological knowledge. Part II. Modelling validation using gravity and magnetic data inversion[END_REF]. Airborne electromagnetic (AEM) data also can constrain the geology in the subsurface: it provides large-scale detailed resistivity mapping of the subsurface [START_REF] D'ozouville | Extensive perched aquifer and structural implications revealed by 3D resistivity mapping in a Galapagos volcano[END_REF][START_REF] Viezzoli | Surface watergroundwater exchange in transitional coastal environments by airborne electromagnetics: the Venice lagoon example[END_REF]. The resistivity variations can be directly used in a 3D geological modelling [START_REF] Jørgensen | 3-D geological modelling of the Egebjerg area, Denmark, based on hydrogeophysical data[END_REF][START_REF] Høyer | 3D geological modelling of a complex buried-valley network delineated from borehole and AEM data[END_REF][START_REF] Sapia | 3D modeling of buried valley geology using airborne electromagnetic data[END_REF]. However, resistivity models deriving from AEM data suffer from uncertainties for depth positioning of geological interfaces, inherent to the sensitivity of the method and equivalence issues. Thus, we propose an alternative method to complement geological data with AEM data conversely. This method consists in jointly interpolating boreholes• geological constraints and information of local slopes derived from the AEM resistivity model using the potential-field method.

The proposed method is illustrated with an AEM survey conducted near Courtenay, France, in February 2009 by SkyTEM ApS. The survey, requested by the Bureau de Recherche Géologique et Minière (French geological survey), covers 11.2 km × 16.5 km (Fig. 1). Being designed for geological and hydrogeological purposes, it provides a semi-detailed electromagnetic (EM) coverage of the studied area. The survey was flown in the north-south direction with 400 m line spacing. The spacing between each EM sounding along flight lines was around 30 m; this represents approximately 12000 EM soundings all over the area.

The studied area is located within the Paris Basin. Its morphology is that of a Cretaceous chalk plateau overlain by about 15 m of weathering clays with flints, together with local thin Quaternary deposits (Fig. 1). The geology is almost layered. However, accurate knowledge of the depth of the clay/chalk interface over the area is missing, both for hydrogeological (i.e., karstic terrain) and industrial (i.e., clay thickness) needs. Courtenay area is well documented with 129 boreholes (Fig. 1) gathered in the Banque de Données du Sous-Sol. These boreholes document the clay/chalk interface; they were drilled for water purposes or as part of seismic studies (Courtenay area is also an oil field) to characterise the low-velocity zone.

The aim of our work was to propose and evaluate a novel method combining geological information from boreholes with the resistivity model from an AEM survey in order to enhance the geological modelling. First, we interpreted, in each AEM sounding, the resistivity contrast corresponding to the clay/chalk interface to derive a geophysical model of this geological interface. In parallel, we modelled this geological interface by traditional interpolation of the borehole data. Then, we achieved a joint interpolation combining the borehole data and trends deriving from the resistivity model. For this purpose, we sampled the slopes of the geophysical model that we previously derived, and we combined these slopes with the borehole data using the potential-field method. The accuracy of this joint model is compared with that of the pure geological and pure geophysical models for different sets of boreholes. Statistical evaluation of the new joint geophysical-geological model is discussed.

H E L I B O R N E T I M E -D O M A I N E L E C T R O M A G N E T I C D A T A

The SkyTEM is a helicopter-borne time-domain electromagnetic system [START_REF] Sørensen | SkyTEM -a new high-resolution helicopter transient electromagnetic system[END_REF][START_REF] Auken | An integrated processing scheme for highresolution airborne electromagnetic surveys, the SkyTEM system[END_REF] initially developed for hydrogeophysical and environmental investigations by the HydroGeophysics Group at Aarhus University, Denmark. SkyTEM system operates in a dualtransmitter mode. The low moment provides early time data for shallow imaging, whereas the high moment allows measuring later time data for deeper penetration. In Courtenay survey, the low moment had a magnetic moment of approximately 3760 Am² with time gates from about 11 to 115 μs and the high moment had a magnetic moment of approximately 140000 Am² with time gates from 73 μs to 9 ms. For a given magnetic moment and noise level, the depth of investigation (DOI) varies depending on the subsurface conductivity and the bandwidth used [START_REF] Spies | Depth of investigation in electromagnetic sounding methods[END_REF]). In the context of our study, the maximal DOI is about 150 m. This allowed investigating the top layers of the chalk underneath a clayish conductive cover (<50 Ωm; [START_REF] Reninger | Geological environment of karst within chalk using airborne time domain electromagnetic data crossinterpreted with boreholes[END_REF]. Navigation data as GPS position, altitude, and tilts of the transmitter loop were processed using filters described by [START_REF] Auken | An integrated processing scheme for highresolution airborne electromagnetic surveys, the SkyTEM system[END_REF]. As part of an environmental study in an anthropised area, particular attention was paid to properly remove noise from the data. Airborne electromagnetic (AEM) data were processed with a singular value decomposition filter [START_REF] Reninger | Singular value decomposition as a denoising tool for airborne time domain electromagnetic data[END_REF], refined through a manual check, in order to remove natural and cultural noises. This was needed to ensure optimal quality data, as small variations in the response were expected [START_REF] Sørensen | SkyTEM -a new high-resolution helicopter transient electromagnetic system[END_REF]. Data were then inverted using the laterally constrained inversion (LCI) algorithm (Auken and Christiansen 2004;[START_REF] Auken | Piece-wise 1D laterally constrained inversion of resistivity data[END_REF]. The LCI is a pseudo-2D inversion scheme; vertical and lateral (along the flight line) constraints are applied on 1D earth models divided into n layers, each being defined by a thickness and a resistivity. In addition, the altitude of the transmitter is also inverted. Results were first obtained with a smooth inversion (19 layers for each 1D model) in order to pre-evaluate the resistivity distribution and adapt the number of layers to the study area. In the following, the study is based on a five-layer inversion, which better discriminates geological interfaces and resistivity contrasts than a smooth inversion (Fig. 2). The DOI was also determined as part of the five-layer inversion (Christiansen and Auken 2012). DOI evaluates the reliable range of depths for each 1D inverted resistivity model. Depths exceeding the DOI were therefore removed from the AEM inversion results.

I M P L E M E N T E D M E T H O D O L O G Y

The clay/chalk interface is well contrasted in the resistivity model (clays are c.a. 20 m in average and chalk is c.a. 100 m in average; see Fig. 2). We were consequently able to interpret an appropriate resistivity contrast corresponding to the geophysical clay-chalk interface comparing the airborne electromagnetic (AEM) resistivity model and boreholes [START_REF] Reninger | Geological environment of karst within chalk using airborne time domain electromagnetic data crossinterpreted with boreholes[END_REF]. The altitude of this contrast was then extracted from each electromagnetic (EM) sounding and mapped throughout the area using a minimum curvature interpolation (Figs. 3 and4). This surface, corresponding to a geophysical approximation of the geological interface, is called "geophysical model" in the remainder of the text.

Since the geophysical model derives from an inversion, it suffers from uncertainties and possible drifts. However, we postulate that its local variations realistically feature the actual variations of depth of the clay/chalk interface. We therefore investigate an original joint processing of the boreholes and EM information in a single combined 3D geological interpolation: in this method, slopes of the geophysical model serve as trends guiding the interpolation of the geological interface constrained at the boreholes. This interpolation uses the 3D implicit potential-field method of [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation[END_REF]; its geostatistical fundamentals are further developed by [START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method[END_REF]. This modelling method was originally developed to interpolate, in the 3D space, borehole information together with structural orientation data measured in the field to derive 3D geological interfaces. In our approach, instead of field structural data, the slopes of the geophysical model are used as orientation data as input for the modelling.

Slopes of the geophysical model were calculated according to [START_REF] Burrough | Principles of Geographical Information Systems[END_REF], with each slope being defined by its orientation and dip. The computation time of the geological modelling is highly affected by the number of input data. We thus sampled the slopes of the geophysical model in order to reproduce its geometry as accurately as possible but with limited number of slopes. The sampling searched for intervals of distance in which the orientations and dips do not vary significantly and calculated average values within each interval. The geophysical model was described with 2391 slopes, corresponding to a computation time of the geological model of about 40 minutes on a standard PC. Figure 3 displays the sampled slopes superimposed on the geophysical model. The main trends are well sampled; tests were made with more slopes without significant changes in the results presented in the following. Slopes were then introduced, as orientation data, in the geological modelling: we computed a combined geophysical-geological geometry of the clay/chalk interface by jointly interpolating the geological contact documented in boreholes with the geophysical slopes using the implicit 3D potential-field method introduced earlier.

E V A L U A T I O N O F T H E J O I N T M O D E L L I N G

In order to evaluate the performance of the proposed joint modelling, we compared three different ways to model the clay/chalk interface over the study area:

r a pure geophysical model, defined by a resistivity contrast that matches the clay/chalk limit documented in boreholes; The accuracy of each of the three obtained models was evaluated as follows. We considered several sets of boreholes as input for the geological and geophysical-geological modelling. Each set contained all 129 boreholes referenced in the area, randomly split into two subsets: constraining and validation boreholes. Constraining boreholes are the boreholes taken into account in the modelling, and validation boreholes are those used to evaluate the models (and therefore left out during the modelling). Then, we computed a geological and a joint model for each set of boreholes. It must be noted that the geophysical model, which we define as a resistivity contrast, remains the same throughout the experiment. For each set of boreholes, the geophysical model and the associated geological and joint models were then compared with the actual clay/chalk interface documented at the validation and constraining boreholes. For each model, the maximum and minimum of the absolute values of the observed errors between the modelled and the actual clay/chalk interface and their average, median, and standard deviation were considered as accuracy indicators. However, preliminary analysis revealed that these indicators display almost the same trends; hence, for simplicity, we present only the absolute average error of each model. Since the number of boreholes is a key parameter in traditional geological modelling, we first considered a test with 10% of validation and 90% of constraining boreholes. Five different sets of boreholes were considered.

Figure 4 presents the three models for one set of boreholes. Obviously, the geological model (Fig. 4b) is the smoothest. It only maps large wavelength variations of the clay/chalk interface, except where the interface intersects the topography, mainly along the two E-W valleys (to the north and south of the study area): along these valleys, the resolution of the geological model directly derives from the resolution of the digital elevation model used. In comparison, the joint model (Fig. 4c) is much more detailed, looking like the pure geophysical model (Fig. 4a) but with less very local variations. The average errors for each of the three models at both the validation and the constraining boreholes for the five sets are plotted in Fig. 5.

At the constraining boreholes (Fig. 5b), the accuracy indicators provide information on the intrinsic accuracy of the models. On the one hand, for the geophysical model, this highlights its average accuracy over the study area; on the other hand, for the geological and joint models, this depends directly on the considered nugget effect when interpolating. For the geological model, the nugget accounts for the uncertainty in the borehole description: depths of geological interfaces documented in boreholes suffer from several uncertainties related to the drilling process and operator interpretation. Moreover, borehole data reflect very local information, and the clay/chalk interface can vary abruptly under the effect of dissolution processes (e.g., [START_REF] Quesnel | Cartographie numérique en géologie de surface -Application aux altérites à silex de l'ouest du Bassin de Paris[END_REF]. For the joint model, the nugget accounts for the uncertainty on the geophysical slopes in addition to the uncertainty in the borehole description. In Fig. 5b, the pure geological (crosses) and geophysical (squares) models are unsurprisingly the most and the less accurate, respectively. Although constrained at the boreholes, the joint modelling (dots) came out more accurate than the geophysical model but does not reach the level of accuracy of the geological model. As the two constraining data (boreholes and geophysical slopes) are of different nature and provide information at different scales, they do not perfectly match and therefore were given a slightly higher nugget.

At the validation boreholes (Fig. 5a), the accuracy indicators provide information on the predictivity of the models. In Fig. 5a, we observe large variations of the average error from one set of boreholes to another since this indicator is calculated on 13 data only (i.e., 10% of the total). The predictivity of the pure geophysical model (squares) is better than that of the pure geological one (crosses). Besides being the least accurate at validation boreholes, the accuracy of the geological model appears to vary greatly over the area (i.e., it is highly dependent on the considered dataset) despite the use of a substantial number of boreholes (i.e., 116). Taking into account geophysical results in the modelling in order to model the clay/chalk interface over the area already proves to be justified at that early stage. Thus, the predictivity of the joint model (dots) is better than that of the geological model and varies less over the area. It is, on average, close to that of the geophysical model and better in some cases (the first and third sets of boreholes in Fig. 5a). As the average error is calculated at 13 validation boreholes only, this result highlights the ability of the joint modelling to locally improve the accuracy of the modelled clay/chalk interface, as compared with the geophysical data.

This first test shows that the joint modelling appears to be the best solution to model the geological interface. Indeed, the proposed method combines the advantages of both the pure geophysical and geological models. We illustrate in Fig. 6a the behaviour of each interpolation along one profile. Boreholes located within a distance of 100 m of the profile are also displayed (constraining boreholes are marked by a C). On the one hand, the geological model (green curve) has a completely different geometry from the geophysical (blue curve) and joint (red curve) models and does not properly predict the altitude of the clay/chalk interface at validation boreholes (without C). On the other hand, the geophysical-geological model, through the adjustment of the geophysical model on constraining boreholes, looks qualitatively like the geophysical model but is closer to the clay/chalk interface documented in boreholes.

These results show the contribution of the joint modelling in comparison with common geophysical and geological models. However, there are two limitations to the previous analysis. First, it was achieved with more than 100 constraining boreholes, which is substantial in comparison with others geo-environmental studies, and second, the model predictivity was evaluated with 13 validation boreholes only. We thus investigated the influence of decreasing the number of constraining boreholes: 90%, 70%, 50%, 30%, and 10% of constraining boreholes were considered successively. Accordingly, models were evaluated with increasingly more validation boreholes. Figure 7 displays the calculated average error of each of the three models at the validation boreholes for the five tested percentages of boreholes. In order to reduce possible low sampling instabilities, for 90% and 70%, the estimator was averaged from five and two sets of boreholes, respectively. At validation boreholes, the average error is the lowest for the joint modelling whatever the number and location of constraining boreholes (Fig. 7). Figure 8 displays the clay/chalk interface mapped by the three models for 10% of constraining boreholes.

As expected, with few constraining boreholes, the geological model comes out smoother (Fig. 8b) and with a large error in comparison with the geometry of the clay/chalk interface as documented in boreholes (Fig. 7). In contrast, the joint modelling keeps a good level of detail, close to the geophysical model. Even with less constraining boreholes, the geometry of the joint model (see the red curve in Fig. 6b) remains closer to the actual clay/chalk interface documented in validation boreholes, though less accurate than when using more constraining boreholes (Fig. 6a). The joint modelling is therefore much less dependent on the number and the location of boreholes taken into account in the modelling than the pure geological modelling. However, the more constraining boreholes we have, the more accurately we can locally adjust the model.

Figure 9 presents the altitude difference between two joint models computed with 90% and 10% of constraining boreholes. The difference between the two models, while small in average (i.e., between +2 and -2 m), can locally increase, where constraining boreholes have been taken into account only in one of the two models.

We computed models from five new sets of boreholes with 10% of constraining boreholes (i.e., 90% of validation boreholes). As the joint modelling is slightly dependent on the number of constraining boreholes (the model being still reliable with few boreholes), this allowed evaluating the accuracy of the joint model at more validation boreholes and so getting a better evaluation of the predictivity of the joint model over the study area for various borehole locations. The calculated average errors of each of the three models at the validation boreholes for the five sets of boreholes are plotted in Fig. 10.

At the validation boreholes, we see that the joint model is close to the geophysical model but is again the most accurate model whatever the location of validation and constraining boreholes considered. This accuracy improvement in comparison to the geophysical model can be explained by the ability of the joint modelling to locally better predict the altitude of the clay/chalk interface (Fig. 5), the joint modelling being efficiently adjusted by only few boreholes. Based on these results, we calculated the median of the error absolute values, in percentage, of each joint model at the 90% validation boreholes. The joint modelling allows predicting the clay/chalk interface over the entire area (11.2 km x 16.5 km) with an average error of 25% (i.e., at 15 m depth, we have an error of 3.75 m). This result must be interpreted with regard to both the size of the area and the noise inherent to the use of borehole data (drilling processes, operator's interpretation, and interface roughness). Moreover, this result has been obtained with 10% of constraining boreholes only (i.e., 13 boreholes only); we expect that the more boreholes we take into account in the modelling, the more the joint model is accurate (Figs. 6 and9).

Then, we looked at the errors as a function of depth of the clay/chalk interface. Figure 11 displays variations of the observed errors with depth at validation boreholes for both the joint and geophysical models for a set of boreholes with 90% of validation boreholes. This result was obtained calculating the median of the absolute values of the errors contained in an interval of 10 m moving every 5 m.

With only 10% of constraining boreholes, the error of the geological model is large (Fig. 10); therefore, we do not display it in Fig. 11. For shallow depths, the accuracy of the joint and geophysical models decreases. Several causes might explain it. For instance, in addition to the resolution limit of the EM method, borehole description is generally less accurate at shallow depths. Moreover, sharp variations of the interface between surveyed lines cannot be predicted by the interpolation of the clay/chalk interface from the 1D layered resistivity models. At intermediate depths, we have the best prediction of the clay/chalk interface with still a slight improvement of the accuracy for the joint modelling in comparison with the geophysical modelling. Note that the clay/chalk interface is documented at a depth higher than 5 m and lower than 20 m in 64% of the validation boreholes. For greater depths, the accuracy of the geophysical model decreases. This loss of resolution with depth is inherent to the EM method. The joint modelling seems to balance this effect: it tends to decrease uncertainties and possible drifts of the geophysical model.

D I S C U S S I O N A N D C O N C L U S I O N

The proposed geological joint modelling utilises regional-scale subsurface investigation by high-resolution airborne timedomain electromagnetic (EM) surveys. It results in a detailed mapping of geological interfaces. In all tested configurations, it improves what is obtained by pure geological or geophysical modelling. Introduction of geophysical slopes between the boreholes locally controls the interface and allows improved prediction of the clay/chalk interface over the study area. It therefore combines the advantages of geological and geophysical modelling. It also tends to decrease uncertainties and possible drifts inherent to the geophysical results. All these aspects make the method very attractive in environmental studies where highly accurate geological models are needed. Still, some specific technical aspects and perspectives are discussed hereafter.

We used a standard minimum curvature method to interpolate the clay/chalk interface interpreted at each 1D layered resistivity model in order to test the direct contribution of slopes in the geological modelling. A more sophisticated interpolation technique such as kriging was also tested but did not significantly modify the results. However, the kriging technique might be more efficient in the case of a more complicated geology.

Moreover, 2391 slopes were used in this study in order to keep a reasonable computation time (40 minutes). Tests were also performed with 3892 slopes for a computation time of 5 hours. We observed no significant improvement in predictivity and higher interferences between slopes and borehole data at a constraining point. Therefore, it seems important to achieve an appropriate sampling of the geophysical slopes, keeping as few slopes as possible.

We showed that the joint modelling is much less dependent on the number of constraining boreholes than the geological modelling. The method could then be particularly adapted in areas with few boreholes. In a geological context such as Courtenay, the resistivity is well contrasted [START_REF] Reninger | Geological environment of karst within chalk using airborne time domain electromagnetic data crossinterpreted with boreholes[END_REF]) and the geophysical model is rather readily interpretable. In more complicated contexts, where it might be difficult to interpret resistivity contrasts, we are confident that it would be possible to generalise this method extracting ad hoc slopes from the resistivity model to calculate the joint model. This perspective of generalisation of the method still has to be evaluated in detail.

Finally, our modelling method is also an alternative to the joint inversion of EM and borehole data. In this alternative approach, the use of geological information is more delicate (can introduce instabilities) and is hardly achievable when boreholes are far from EM survey lines. Moreover, the effect of constraining boreholes is generally very local. For the purpose of deriving directly geological interfaces from EM surveys, our integrated approach is easier as all boreholes can be integrated and associated together with a standard robust geophysical inversion; the boreholes optimise the modelled geological interfaces locally and on intermediate to long ranges.

A C K N O W L E D G E M E N T S

This work was carried out as part of a PhD co-financed by the Région Centre and the BRGM. Geophysical data on which the study is based were funded by the BRGM and ANR through the Carnot-BRGM label. The authors would like to thank two anonymous reviewers who suggested reworking of several important points of the original manuscript.

Figure 1

 1 Figure1Simplified geological map of the study area (modified from[START_REF] Pomerol | Carte géologique de la France (1/50 000[END_REF]. Flight lines are displayed in gray and boreholes are indicated by black triangles. The black rectangle in the insert map locates the study area, in France.

Figure 2

 2 Figure 2 5 layer resistivity soundings along a flight line; a fading was applied on results below the depth of investigation. Boreholes located within a distance of 50 m of the profile are also displayed: Quaternary deposits in gray, clays in white and chalk in black. Vertical exaggeration is 2.

Figure 3

 3 Figure 3 Sampling of the slopes of the geophysical top of the chalk. Black points and black lines, with lengths proportional to the dip, represent respectively the location and the dip/orientation of the sampled slopes.

Figure 4

 4 Figure 4 Modelling of the clay/chalk interface over Courtenay area with a set of boreholes with 90 % of constraining boreholes; (a) geophysical model, (b) geological model, (c) geophysical-geological model. Constraining boreholes are indicated by black dots.

r

  a pure geological model, obtained by interpolation of the clay/chalk interface in the boreholes, as well as the outcrops of this contact mainly along the two E-W valleys, to the north and south of the study area; r a joint geophysical-geological model, obtained by the joint interpolation of the clay/chalk interface in the boreholes, together with the slopes of the geophysical model, using the implicit 3D potential-field method.

Figure 5

 5 Figure 5 Average errors of the geological (crosses), geophysical (squares) and joint (dots) models for the five sets of boreholes considered; (a) at the validation boreholes, (b) at the constraining boreholes.

Figure 6

 6 Figure 6 Cross-section displaying the geological (in green), geophysical (in blue) and joint (in red) models for (a) 3 and (b) 1 constraining boreholes within the profile. Constraining boreholes are indicated by a C. Boreholes within a distance of 100 m of the profile are displayed; clays are represented in white in boreholes and chalk in black. Short black lines and black dots on the red model display respectively the geophysical slopes oriented along and sub-perpendicular to the profile within a distance of 100 m.

Figure 7

 7 Figure 7 Average errors of the geological (crosses), geophysical (squares) and joint (dots) models at the validation boreholes in function of the percentage of constraining boreholes.

Figure 8

 8 Figure 8 Modelling of the clay/chalk interface over Courtenay area with 10 % constraining boreholes; (a) geophysical model, (b) geological model, (c) geophysical-geological model. Constraining boreholes are indicated by black dots.

Figure 9

 9 Figure 9Difference in altitude between two joint models computed with 90 and 10 % of constraining boreholes. Constraining boreholes location appears as black dots and circles for respectively 90 and 10 %.

Figure 10

 10 Figure 10For 10 % of constraining boreholes: average errors of the geological (crosses), geophysical (squares) and joint (dots) models at the validation boreholes for the five sets of boreholes considered.

Figure 11

 11 Figure 11Evolution of the median of the error with depth of the clay/chalk interface documented in the validation boreholes, for the geophysical (squares) and the joint (dots) models.


Sakina Mahboubi
email: sakina.mahboubi@inria.fr

Ph D Thesis

Reza Akbarinia

Patrick Valduriez

Privacy Preserving Query Processing in the Cloud

INTRODUCTION

Nowadays, cloud data outsourcing provides users and companies with powerful capabilities to store and process their data in third-party data centers. However, when a user stores her sensitive data in a public cloud, they becomes vulnerable to several attacks, e.g., from the employees of the cloud provider.

One solution for protecting the user data against attacks is to encrypt the data before sending them to the cloud servers. Then, the challenge is to answer user queries over encrypted data. A naive solution for answering queries is to retrieve the encrypted database from the cloud to the client, decrypt it, and then evaluate the query over plaintext (non encrypted) data. This solution is not practical, in particular for large databases.

In this PhD thesis, we are interested in processing topk queries on encrypted data. These queries have attracted much attention in several areas of information technology such as sensor networks [START_REF] Yang | An efficient top-k query processing framework in mobile sensor networks[END_REF], stream management systems [START_REF] Wang | SKYPE: top-k spatial-keyword publish/subscribe over sliding window[END_REF][START_REF] Shen | Efficiently monitoring top-k pairs over sliding windows[END_REF] and spatial data analysis [START_REF] Chen | Answering why-not spatial keyword top-k queries via keyword adaption[END_REF][START_REF] Choudhury | Batch processing of top-k spatial-textual queries[END_REF]. A top-k query allows the user to specify a number k, and the system returns the k tuples which are most relevant to the query. The relevance degree of tuples to the query is determined by a scoring function.

There have been many different approaches proposed for processing top-k queries over plaintext data. Two of the best known approaches are FA [START_REF] Fagin | Combining fuzzy information from multiple systems[END_REF] and TA [START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF] that work on sorted lists of attribute values. These approaches, particularly TA, can find efficiently the top-k results because of smart strategies for deciding when to stop reading the database. However, TA, FA and all other efficient top-k approaches developed so far assume that the data are in plaintext, and there is no efficient solution capable of evaluating efficiently top-k queries over encrypted databases.

When we think about top-k query processing on encrypted data, the first idea that comes to mind is the utilization of a fully homomorphic encryption cryptosystem, e.g. [START_REF] Gentry | Fully homomorphic encryption using ideal lattices[END_REF], which allows to do arithmetic operations over encrypted data. this type of encryption allows to compute the overall score of data items over encrypted data. However, fully homomorphic encryption methods are very expensive in terms of encryption and decryption time. In addition, they do not allow to compare the encrypted data, and to find the top-k results.

We proposed efficient approaches, called EncFA and Buck-Top, for processing top-k queries over encrypted data. We evaluated their response time over encrypted data with that of the TA algorithm over original (plaintext) data. Our results show that the response time of BuckTop over encrypted data is close to TA over plaintext data, and even better over some large databases.

In the rest of this paper, we fist define the problem we address. Then, we briefly introduce our proposed approach, and then we discuss the related work.

PROBLEM DEFINITION

The problem which we address is top-k query processing over encrypted data. Let D be a database, and e(D) be its encrypted version such that each data c ∈ e(D) is the ciphertext of a data d ∈ D, i.e. c = Enc(d) where Enc() is an encryption function. The database e(D) is stored in a remote server.

Given a number k and a scoring function f , our goal is to develop an algorithm A, such that when A is executed over the database e(D), its output contains the ciphertexts of the top-k results, i.e. those that can be found by executing a correct top-k algorithm over the database D.

APPROACH

The architecture of our system for query processing over encrypted data is composed of two parts: Trusted client and Service provider.

Trusted client is responsible of user data encryption and decryption, user access control and secure key management. When a query is issued by a user, the trusted client checks the access rights of the user. If the user does not have the required rights to see the query results, then her demand is rejected. Otherwise, the issued query is transformed to a query that can be executed over the encrypted data and is sent to the service provider. the Trusted client decrypts the received results (calculated by the service provider) and returns to the user the k data item which are the response of the query launched.

Service provider stores the encrypted data, executes on them the algorithms provided by the trusted client, and returns the results to it.

We propose two approaches for processing top-k queries over encrypted data. The first approach, called EncFA, uses a probabilistic (semantic) encryption scheme to encrypt the attribute values of the database items. It also uses a deterministic scheme for encrypting the IDs of data items. Then, it sends the encrypted database to the service provider. the data items after encryption must have the same order that they had before being encrypted.

The second approach, called BuckTop, is more efficient than EncFA. It uses a bucketization technique to partition the data items in each list into a set of buckets. also it uses tow types of encryption schema to encrypt the database; one is deterministic used to encrypt data item IDs and the other is probabilistic used to enrypt the attribute values of the data items. BuckTop includes a top-k query processing algorithm that works on the encrypted data of the buckets, and returns a set containing the top-k results. It also includes an efficient filtering algorithm that filters the false positives as much as possible in the server. We prove theoretically the correctness of the BuckTop approach.

RELATED WORK

A first important paper in top-k query processing is [START_REF] Fagin | Combining fuzzy information from multiple systems[END_REF], which models the general problem of answering top-k queries using databases organized into lists of data items sorted by their local scores. One of the most efficient algorithms over sorted lists is the TA algorithm, which was proposed by several groups [START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF][START_REF] Güntzer | Towards efficient multi-feature queries in heterogeneous environments[END_REF]. However, all these algorithms assume that the data scores are available as plaintext, and not encrypted.

There have been also some proposed solutions for secure kNN query processing, e.g. [START_REF] Elmehdwi | Secure k-nearest neighbor query over encrypted data in outsourced environments[END_REF][START_REF] Choi | Secure knn query processing in untrusted cloud environments[END_REF][START_REF] Yao | Secure nearest neighbor revisited[END_REF]. The problem is to find k points in the space that are the nearest to a given point. This problem should not be confused with the top-k problem in which the given scoring function plays an important role, such that on the same database and with the same k, if the user changes the scoring function, then the output may change. Thus, the proposed solutions proposed for kNN cannot deal with the top-k problem.

The bucketization technique has been used in the literature for answering range queries over encrypted data, e.g. [START_REF] Hore | Secure multidimensional range queries over outsourced data[END_REF][START_REF] Li | A data-and workload-aware query answering algorithm for range queries under differential privacy[END_REF] where Hore et al. use this technique, and propose optimal solutions for distributing the encrypted data of a database to the buckets in order to guarantee a good performance by reducing the number of false positives while preserving a high security level. The techniques developed in [START_REF] Hore | Secure multidimensional range queries over outsourced data[END_REF][START_REF] Li | A data-and workload-aware query answering algorithm for range queries under differential privacy[END_REF] can be used in our system for an optimal distribution of the encrypted data in the buckets.

The only paper which we found about top-k query processing over encrypted data is [START_REF] Xianrui Meng | Declarative cleaning of inconsistencies in information extraction[END_REF] published in arXiv.org. The proposed architecture assumes the existence of two noncolluding servers s1 and s2 in two different clouds. One of the servers, say s2, has the decryption keys, and the other one, say s1, stores the data. Top-k query processing proceeds by using the TA algorithm and accessing the encrypted data in s1, such that after reading each data in s1, its encrypted local scores are sent to the server s2 (using a special protocol) where they are decrypted and compared with the TA threshold. Our assumptions about the cloud are different. In our solution, we do not need to trust on any remote server, and during the top-k query processing, we do not decrypt the encrypted data in the cloud servers. In addition, the solution in [START_REF] Xianrui Meng | Declarative cleaning of inconsistencies in information extraction[END_REF] needs a lot of communications between remote servers (i.e., at least two messages after each sorted access). This solution that is not efficient and incurs a high latency in the query processing time.

to the best of our knowledge, in the literature there is no efficient solution for processing top-k queries over encrypted data. In this work, we propose such a solution.

(c) 2016 ,

 2016 Copyright is with the authors. Published in the Proceedings of the BDA 2016 Conference (15-18 November, 2016, Poitiers, France). Distribution of this paper is permitted under the terms of the Creative Commons license CC-by-nc-nd 4.0. (c) 2016, Droits restant aux auteurs. Publi dans les actes de la confrence BDA 2016 (15 au 18 Novembre 2016, Poitiers, France). Redistribution de cet article autorise selon les termes de la licence Creative Commons CC-by-nc-nd 4.0. BDA 2016, 15 au 18 Novembre, Poitiers, France.