Frequency Domain Electromagnetic mapping for delineating subsurface structures related to the historical port of Emporiae

In recent years, there has been a growing interest in the application of geophysical methods to reconstruct the palaeo‐landscapes of sites of special historical interest in support of the planning of archaeological researches. Given the extent of the surface to be investigated, electromagnetic methods have proven to be very suitable for their speed, resolution and versatility for this objective. In particular, coastal areas of the Mediterranean have undergone significant changes in the position of the coastline, because of changes in sea level and sediment inputs that have covered natural harbours used for the establishing colonies.

different port areas, each with its own specific characteristics and adapted to specific needs.
This multidisciplinary project is structured in two main research lines that are completely complementary. The first refers to the archaeological research and involves the excavation in different sectors of Empúries, as well as the study of all the materials recovered in these surveys. The second deals with the geomorphological and palaeoenvironmental reconstruction of these ancient port areas, since knowledge of this is essential to contextualize the data obtained in archaeological interventions. The deep transformation of the landscape around Empúries requires historical interpretation to keep in mind the constant interrelationship between the dynamics of occupation and the dynamics of the landscape.
The study of this historical period focuses most of the interventions and other actions contemplated in this project. The origin and evolution of Emporion is therefore presented as a result of several factors, among which the availability of various port spaces is of particular significance. The case of the so-called "natural harbour" is one of these cases because it was perfectly adapted to the kind of trade developed in the Emporion. Besides, it allowed direct communication with inland territories due to its connection with the estuary of the Fluvià and the northern arm of the Ter.
The archaeological activity during 2018 and 2019 focused on the excavation of the northern quarter of the Greek city, which bordered the natural harbour. This work has recovered the original topography, revealing the cliff boundary that limits with the small harbour bay, above which was the Greek settlement. The original topography had been greatly altered because of the intense sand filling process, which during the late ancient period completely refilled the natural harbour.
Communication between the city and the port was made through a ramp located at the northwest corner, where the rock had a softer relief and where there was a small beach. The urban layout of the Greek centre perfectly adapted to the natural topography of the land, which was markedly higher than the port and the coastline.
These archaeological interventions have been complemented with new geophysical surveys in the interior of the natural harbour, as well as various geological surveys, with the aim of gaining a deeper understanding of the characteristics of this port space and establishing the chronology and causes of the subsequent filling process.
Previous researches about sedimentology and evolution of coastline suggested that alluvial sediments of the Fluvià River and sands transported by the sea concealed the historical port of Emporiae (Marquès & Julià, 1983).
The current name of Empúries comes from the Greek term Emporion, which means shopping centre or mall, and faithfully described the purpose of the site, because the city was built initially in the delta of the Fluvià River, crossing several trade routes, and with a natural port, which offered adequate protection to commercial ships.
The first Greek settlement is dated to the 6th century BC. The location was in the present town of St. Martin Empúries (Palaiapolis) which at that time was surrounded by water, and later was integrated to the coast by an isthmus thanks to the contributions of the sediments of the Fluvià River. In the 5th century BC the Greeks moved Emporion, from St. Martí Empúries to the south, at an emerged hill in the shoreline (Neapolis). Then, Emporion quickly became one of the most important commercial ports in the Mediterranean ( Figure 1).

| HISTORICAL FRAMEWORK
The various nuclei of population that, over time, were found in Empúries and its more immediate environment show the importance that this commercial place had during ancient history. During the 7th century BC, the ideal location of the Empordà coast favoured the stable occupation of the various hills in the Empúries area by indigenous communities and the establishing of initial commercial contacts, which soon after, were to culminate in the setting up of the well-known Greek Phocaean Emporion on the promontory of Sant Martí d'Empúries, the Palaia Polis.
The particular morphology and topography of the coastal area, with a small natural bay that served as a port between the island/ peninsula of Sant Martí to the north, and the much more wider coastal promontory to the south, allowed that, shortly after the first settlement (Paleapolis), a new Greek city (Neapolis) was built, which kept the name of Emporion.
The origin and evolution of Emporion were always linked to its commercial vocation and its port, as can be seen by the fact that centuries later it also became the Roman gateway to the Iberian Peninsula, initially to solve the war conflict with the Carthaginians and later to contribute to the control and conquest of Hispania. The strategic importance of Empúries at that time is evident from the setting up of a permanent Roman military camp under the shelter of the Greek nucleus, which, once abandoned, was used as a base to create a new Roman city in the early 1st century BC. The settlement was orthogonal in shape and covered a surface area of 22.5 ha, most of which have yet to be discovered. Towards the change of the era, the Greek and the Roman nuclei merged into a single site that we know as municipium Emporiae ( Figure 2).

| METHODOLOGY
Field geophysical surveys started in 1996 in the framework of a project leaded by the German Archaeological Institute using a set of different geophysical methods: Vertical electrical soundings (VES), ground probing radar (GPR) and magnetics (Blech et al., 1998). VES showed roughly the structure of the basin and were complemented by a borehole, while GPR and magnetics did not provide results of interest. Subsequently, an ERT survey was carried out but we face logistical difficulties in accessing because it is an agricultural area and the slowness of data acquisition. The location of the previous geophysical measurements is shown in Figure 6. Recently our efforts has been concentrated on frequency-domain electromagnetic method (FDEM) that has shown to be very effective for similar archaeological researches (Bates & Bates, 2000;Tabbagh, 1986).
According to the objectives of the present study, the geophysical surveys were carried out through gradual stages. First, VES for 1D modelling, ERT for 2D modelling and finally FDEM for mapping.

| FDEM
The assessment of soil electrical conductivity properties using low frequency induction methods (Keller & Frischknecht, 1966;Wait, 1955) has decades of many successful applications for environmental (Nobes, 1996), hydrogeological (Boaga, 2017) and archaeological studies (Christiansen et al., 2016) among others. In this method, a primary electromagnetic field is induced by the transmitting coil carrying a time-varying electric current at a set frequency, which generates a (primary) magnetic field into the subsurface. The resulting secondary magnetic field is measured together with the primary magnetic field at the receiver coil. The ratio between the secondary and primary magnetic field is recorded as in-phase and quadrature-phase data ( Figure 3). The apparent electrical conductivity (ECa) of the bulk soil can be obtained through the quadrature-phase as a depth weighted electrical conductivity value of a homogeneous medium using the formula given by McNeill (1980): where f is the frequency (Hz), s is the coil separation, μ o is the magnetic permeability of free space (4π10 À7 H/m) and (H s /H p ) Qu is the Q u component of the secondary H s to primary H p magnetic field coupling ratio (McNeill, 1980). The formula is an approximation based on the assumption of operating the instrument at low induction number (LIN) and null clearance above the ground conditions. The dimensionless LIN parameter is defined as the ratio between the instrument coil separation and the skin depth δ, where the skin-depth in turn is defined as the distance wherein a plane wave is attenuated by 1/e (approx. 37%) of the value at the surface (Spies, 1989).
A Geonics EM31 MK2 ground conductivity meter has been used.
The background conductivity noises of the device are 0.1 mS/m for Data were directly stored to a DL600 data-logger. Repeated measurements were recorded at the beginning and end of each profile to perform corrections for intersection errors (miss ties) that could affect the homogeneity of data during gridding.

| VES
DC resistivity surveys in the study area were carried out in the framework of project held in 1996 by German Archaeological Institute project (Blech et al., 1998). Five vertical electrical sounding (VES) were recorded along a profile crossing the study area ( Figure 6). This preliminary survey was conducted to define thickness of the sediments over the high resistivity basement formed by cretaceous limestones.
The geophysical equipment used was a DataRes digital resistivity meter from Ambrogeo. A Schlumberger array was selected for this survey because his suitable sensitivity for mapping shallow electrical resistivity variations. The distances between current electrodes (AB) were ranged logarithmically from 1 m to 60 m. This spread was considered long enough to penetrate the sedimentary cover to a depth of about 15 m.
The electrode coupling with the ground was checked before carrying out the measurement for each of the four electrodes at F I G U R E 2 Location of the different historical sites referred in the text [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 3 Schematic representation of FDEM principles. H p is the primary field generated the transmitter coil T x ; Hs is the secondary field recorded at the receiver coil Rx; dotted lines are Eddy currents (modified from Boaga, 2017) [Colour figure can be viewed at wileyonlinelibrary.com] each measurement. The IX1D software (Interpex, 2008) was used to process the VES data sets and obtain the subsurface resistivity model. The inversion program incorporates a ridge regression routine to achieve a best fit to the observed data (Inman, 1975). The best fit can be achieved either automatically through the software or manually by adjusting the layer thickness and resistivity values.
To overcome the ambiguity problem in the inversion of resistivity curves, the lithologic information of a borehole was used to build the input models of the inversion process. The borehole S1 was drilled close to VES2 (Figure 6), and some model parameters were fixed during the inversion process according to the available geological data. Apparent resistivity values were first checked for erroneous values that may occur mainly due to bad electrode coupling with the ground that show up as spikes in the apparent resistivity. For all lines, we did then 2D inversion using ram Res2DInv ® (Geotomo Software, 2010) including the topography. During the inversion process, the root-mean-square value of the difference between experimental data and the updated model response is used as a criterion to assess the convergence (Loke & Barker, 1996). For most lines a root mean square (RMS) data adjustment better than 3% was obtained.

| FDEM
Apparent electrical conductivity maps have been obtained from FDEM values after statistical data processing and gridding using Krigging interpolation method and taking into account the anisotropy of sampling (Hansen, 1993). The result was respectively a 5 Â 5 m grid for any coil orientation. Then, colour contour maps have been

| VES
Experimental apparent resistivity data from all VES curves clearly showed three-layer H type models (Koefoed, 1979) with a thin high resistivity top layer, interpreted as a soil cover, an intermediate layer of fine size sand mixed with clay of low electrical resistivity and thicknesses ranging from 5 to 10 m and a high resistivity layer at the bottom interpreted as the Cretaceous limestone. The inverted model of VES 2, constrained by the geological log from the borehole S 1, drilled few meters apart, is shown in Figure 8. The thickness of the first layer kept fix to 1.5 m, and the thickness of the second layer was considered the free variable.
Taking into account that some authors (Callegary et al., 2007;Christiansen et al., 2016;Deidda et al., 2014;Reid & Howlett, 2001, among others)  The response of the model was calculated following the formula given by McNeill (1980): where σ a is the bulk apparent conductivity of the layered earth, σ 1 , σ 2 and σ 3 are the electrical conductivities of the top, intermediate and bottom later, z 1 and z 2 the depths to the second layer and third layer, respectively. R V or R H are the relative contributions of vertical and horizontal components to the secondary magnetic field or apparent conductivity from all material below a depth z given by Cross-section of the ERT profile which location is showed on Figure 6 [Colour figure can be viewed at wileyonlinelibrary.com] The functions ;(z) and R(z) define the relative influence of current flow as a function of depth a z is the depth divided by the intercoil space. Always following McNeill (1980), the functions R V (z) and R H (z) can be easily calculated by Then, FDEM apparent electrical conductivity data close to borehole S1 have been inverted keeping invariant all the parameters of the model except the thickness of the second layer (Guérin et al., 1996).
The apparent electrical conductivity of the model for vertical dipoles was 55.80 mS/m whereas the observed apparent conductivity was 56.2 mS/m.
The same 1D inversion process has been applied to all FDEM conductivity values and an iso-depth map of the basement has been generated ( Figure 10). The iso-depth map shows the existence of structures that can be assumed as platforms and channels. This interpretation is following the hypothesis advanced by the archaeologists that is reproduced ( Figure 11) and with a borehole drilled that found the limestone substrate at 5.20 m depth. Besides, the radiocarbon dating of plant remains at 4.80 m depth gave an age of 2020 ± 70 years, which is consistent with the historical period of the port. This area, which is currently outside the archaeological complex, in ancient times, was part of the same urban centre. At some points, the stratigraphic sequence made it possible to identify different successive levels of occupation, as well as some sediments that directly covered the natural rock. The ceramic materials associated with these strata confirm that the occupation of this sector occurred during the initial stage of the Neapolis, starting in the second half of the 6th century BC (Marcet & Sanmartí, 1989).

| CONCLUSIONS
The results obtained in this study indicate that shallow electromagnetic induction is a cost-effective alternative for mapping the buried palaeo-landscape related to harbours and coastal plains laying over a high resistivity basement. In our case, the existence of a thin low electrical conductivity layer overlaying a fine sand mixed with clay, the resolution of subsurface parameters is much better than in the opposite case of a highly conductive layer overlaying a moderately conductive layer.
Particularly, the geophysical images and resulting interpretation from data obtained in this study confirmed the existence of a well-defined basin between the Greek archaeological sites named Palaiapolis and Neapolis.
Then, we can conclude that the basin was used in Antiquity as a natural port as has been suggested by the Archaeologist from long ago without conclusive evidence. Nevertheless, further archaeological diggings and geophysical surveys should be conducted.