Multi-Sensor Geomagnetic Prospection : A Case Study from Neolithic Thessaly , Greece

Multi-sensor prospecting is a fast-emerging paradigm in archaeological geophysics. Given suitable ground conditions for navigation, sensor arrays drastically increase efficiency in data collection. In particular, geomagnetic prospecting benefits from this development. Despite these advancements, data processing still lacks a best-practice approach. Conventional processing methods developed for gridded data has been challenged by sensor arrays “roaming” in the landscape. In realization of the issue, the Innovative Geophysical Approaches for the Study of Early Agricultural Villages of Neolithic Thessaly (IGEAN) Project explored various innovative techniques for the betterment of the multi-sensor geomagnetic data processing. As a result, a modular pipeline is produced with minimal user intervention. In addition to standard steps, such as data clipping, various other algorithms have been introduced. This pipeline is tested over 20 Neolithic settlements in Thessaly, Greece, three of which are presented here in detail. The proposed workflow provides drastic improvements over raw data. As a result of these improvements, the IGEAN project revealed astonishing details on architectural elements, settlement enclosures, and paleolandscapes, changing completely the existing perspective of the Neolithic habitation in Thessaly.


Introduction
The Neolithic period in Europe was a key epoch.Migrant groups gave way to permanent agrarian societies occupied with animal husbandry and the cultivation of food-crops for sustenance.Conceptualization of the Neolithic surroundings became viable for the first time in durable architectural forms.Thessaly (Central Greece) is of critical importance in this transformation serving as the gate to what would become the widespread Neolithization of Europe, which irreversibly altered the course of human history [1].The high density of Neolithic activity and mounded settlements (locally called magoules; magoula (sng.)), the early date of occupation, and the multiple trajectories that it followed make Thessaly a crucial region for exploring and comprehending the pathways in which Neolithic culture emerged and progressed in the European continent (Figure 1).This paper focuses on the geomagnetic prospection component of the IGEAN and reveals three prime examples from the project.In doing so, it aims to show the power of multi-sensor systems as an emerging paradigm in archaeology.Given that the surface conditions are appropriate for extensive data collection, the efficiency of multi-sensors over single or dual sensor systems is beyond doubt.It is now truly possible to conduct "geophysical surveys as landscape archaeology" [2].

State-of-the-Art
There are two main handheld gradiometers, deployed for archaeological prospecting (Geoscan FM series and Bartington Grad601/DualGrad601).These sensors should be carefully balanced before surveys in order to acquire the best available geomagnetic data.Due to various reasons (such as temperature fluctuations and abrupt shocks) these sensors usually drift throughout the survey and requires repeated trimming.The operator gait is also an important parameter in the collection of data and obstacles or user errors may require repetition of parts of the survey.Finally, these common sensors require establishing collection grids prior to any survey.The grids are further divided into equally-spaced traverses on which the surveyor walks back and forth at a (relatively) constant speed in order to obtain an even distribution of geomagnetic data.Within a new scientific agenda of investigating Neolithic cultures, the Laboratory of Geophysical-Satellite Remote Sensing and Archaeo-environment (GeoSat ReSeArch Lab) of the Foundation for Research and Technology-Hellas (FORTH) initiated (2013-2015) a remote sensing campaign to study the physical landscape and dynamics of Neolithic settlements within the coastal hinterlands of Eastern Thessaly.The IGEAN (http://igean.ims.forth.gr/)research project is the first systematic and extensive geophysical survey of Neolithic magoules in the region and is the first to utilize multi-sensor geomagnetic sensors in the Eastern Mediterranean.
This paper focuses on the geomagnetic prospection component of the IGEAN and reveals three prime examples from the project.In doing so, it aims to show the power of multi-sensor systems as an emerging paradigm in archaeology.Given that the surface conditions are appropriate for extensive data collection, the efficiency of multi-sensors over single or dual sensor systems is beyond doubt.It is now truly possible to conduct "geophysical surveys as landscape archaeology" [2].

State-of-the-Art
There are two main handheld gradiometers, deployed for archaeological prospecting (Geoscan FM series and Bartington Grad601/DualGrad601).These sensors should be carefully balanced before surveys in order to acquire the best available geomagnetic data.Due to various reasons (such as temperature fluctuations and abrupt shocks) these sensors usually drift throughout the survey and requires repeated trimming.The operator gait is also an important parameter in the collection of data and obstacles or user errors may require repetition of parts of the survey.Finally, these common sensors require establishing collection grids prior to any survey.The grids are further divided into equally-spaced traverses on which the surveyor walks back and forth at a (relatively) constant speed in order to obtain an even distribution of geomagnetic data.
A multi-sensor array is a significant improvement over handheld gradiometers.There is no need for manual trimming of individual sensors; software routine calculates the adjustments needed to the readings to 'balance' the signal.Furthermore, the connection between sensor arrays and accurate navigation systems makes it possible to "roam" on the landscape and collect high resolution data in an efficient manner; i.e., with the advancement of the multi-sensor technology there is no need for establishing survey grids anymore.Nevertheless, such an extensive mapping with an extremely dense point pattern brings about new challenges.Navigated through a Global Navigation Satellite System (GNSS) differential unit, a large volume of data is collected in a non-gridded way, having areas that are overlapped by different paths and by different sensors, creating various types of noise to the raw data.The particular paper addresses ways that deal with these problems and proposes a more automatic way of producing quality magnetic maps.

Study Area and Previous Research
Thessaly is a closed geographical unit with mountainous borders and an eastern coastline with access to the Aegean Sea.Its interior consists of several fertile basins at different altitudes created by tectonic episodes and the hydrological activity of its rivers and lakes.Two major floodplains, the Karditsa Plain to the west and the Larissa Plain to the east, are separated by the Revenia Hills (~700 m).In spite of the orography of the region, communication between the eastern and western plains was common among Neolithic peoples judging from the wide distribution of material culture [3][4][5].Thessaly contains remarkable evidence of continuous habitation in all phases of the Neolithic period, documented mainly in the form of magoules (~350 have been identified) attracting the interest of scientific inquiries since the beginning of the 20th century and after [6][7][8][9], partly due to the outstanding state of preservation and monumental character of some sites (diameter >200 m and height ~8 m).In the 1980s and 1990s, research turned to the identification of new sites (largely through rescue excavations) and to the delineation of Neolithic settlement patterns [10,11].Recently, Alexakis et al. [12,13] focused on the creation of an updated corpus of sites, a reconstruction of the geomorphology of Thessaly, the use of space-borne imagery to identify previously unknown sites in the region and GIS spatial analysis to characterize the dynamics of Neolithic habitation.
Despite the significance of the area for the studies on Neolithic Europe, key questions are still pending to be scrutinized.These include: environmental factors behind the establishment of early farming communities, spatial development of agricultural villages, land use patterns, and farming strategies around sites, dynamics of ecological niches, internal organization of structures within different types of settlements, and the meanings behind settlement enclosures.The above constituted key questions for the IGEAN project as it was the first time that geophysical surveys were able to scan extensively 24 Neolithic settlements and try to identify critical characteristics of their internal organization and of their surroundings (Figure 1).Despite of the manifold approach used in the survey of the Neolithic settlements, it was the multisensory magnetometry that managed to provide the bulk of the information related to the above structural features of the sites.

Instrumentation
As part of the IGEAN project, a SENSYS Sensorik and Systemtechnologie GmbH (Bad Saarow, Germany) MX Compact was used for geomagnetic data collection (Figure 2).The SENSYS MX is a multi-channel system equipped with FGM650/3 single-axis fluxgates with a standard measurement range of ±10,000 nT.The system is dependent on accurate positioning information with the (National Marine Electronics Association) NMEA 0183 standard, synchronized to one point-per-second (PPS).
A proprietary software (MonMX, version 4.0) helps the surveyor in navigation, using GNSS data and predefined sensor array geometry.Gradiometers are mounted on a non-magnetic carrier which can hold up to 16 sensors with 25 cm separation.Depending on the needs of the geophysical survey and ground conditions a numerous carrier-sensor combination is also possible.In the initial setup, the distance between sensors is fixed along the array (x-axis).However, measurements along the track (y-axis) depend on the speed of survey since the sensors register magnetic data in fixed time intervals; a faster movement of array results in fewer readings, and vice versa.Survey conditions or navigation skills may also result in denser or sparser data at the horizontal plane as well.Therefore, the multi-sensor paradigm dictates non-regularized data collection.While non-regularity hinders to use common geophysical data processing methodologies, such as despiking or empirical drift correction, it also opens up new data exploration possibilities.
In surveying Neolithic settlements, the selected configuration was an eight-sensor setup with 50 cm separation in order to increase the prospection coverage.In a couple of cases, this configuration was modified for shorter arrays, especially in places with intervening obstacles (such as olive trees).A Javad Triumph-1 receiver is used for the accurate positioning of the multi-channel system in the real-time kinematic (RTK) mode.

Data Analysis
For data analysis, a set of customized data processing routines was developed in the MATLAB (MathWorks) computing software.The workflow includes standard steps, such as American Standard Code for Information Interchange (ASCII) import and parsing, data clipping, reporting descriptive statistics, etc., for each survey track and for the entirety of the survey coverage.When required, data stripes are removed using Fast Fourier transformation (FFT).
At this point, it should be further noted that striping remains an issue for multi-sensor arrays.Effectiveness of transformation mainly depends on the modifications of the geomagnetic data spectrum.Therefore, the results of the FFT are still user-determined and the problems inherent to the gridded data are also passed over to the multi-sensor array datasets.
In order to provide further processing of the non-gridded dense point pattern geomagnetic data, new algorithms are also introduced specifically for data sparsing and outlier removal (Figure 3).The data processing workflow proved to be efficient and provided drastic improvements over raw data (Figure 4).Gradiometers are mounted on a non-magnetic carrier which can hold up to 16 sensors with 25 cm separation.Depending on the needs of the geophysical survey and ground conditions a numerous carrier-sensor combination is also possible.In the initial setup, the distance between sensors is fixed along the array (x-axis).However, measurements along the track (y-axis) depend on the speed of survey since the sensors register magnetic data in fixed time intervals; a faster movement of array results in fewer readings, and vice versa.Survey conditions or navigation skills may also result in denser or sparser data at the horizontal plane as well.Therefore, the multi-sensor paradigm dictates non-regularized data collection.While non-regularity hinders to use common geophysical data processing methodologies, such as despiking or empirical drift correction, it also opens up new data exploration possibilities.
In surveying Neolithic settlements, the selected configuration was an eight-sensor setup with 50 cm separation in order to increase the prospection coverage.In a couple of cases, this configuration was modified for shorter arrays, especially in places with intervening obstacles (such as olive trees).A Javad Triumph-1 receiver is used for the accurate positioning of the multi-channel system in the real-time kinematic (RTK) mode.

Data Analysis
For data analysis, a set of customized data processing routines was developed in the MATLAB (MathWorks) computing software.The workflow includes standard steps, such as American Standard Code for Information Interchange (ASCII) import and parsing, data clipping, reporting descriptive statistics, etc., for each survey track and for the entirety of the survey coverage.When required, data stripes are removed using Fast Fourier transformation (FFT).
At this point, it should be further noted that striping remains an issue for multi-sensor arrays.Effectiveness of transformation mainly depends on the modifications of the geomagnetic data spectrum.Therefore, the results of the FFT are still user-determined and the problems inherent to the gridded data are also passed over to the multi-sensor array datasets.
In order to provide further processing of the non-gridded dense point pattern geomagnetic data, new algorithms are also introduced specifically for data sparsing and outlier removal (Figure 3).
The data processing workflow proved to be efficient and provided drastic improvements over raw data (Figure 4).

Outlier Removal Based on Space-Phase Method
Due to the large volume of geomagnetic data collected per unit area, it is not practical to employ standard outlier detection and removal/replacement methodologies.In a 20 m × 20 m survey area, an eight-sensor array configuration with 0.5 m separation may contain 150,000 data points collected at an average walking speed.This is an immense dataset in comparison to 1600 geomagnetic data points collected within a 20 m × 20 m survey grid with a 0.5 m transect and sampling spacing (Figure 5).

Outlier Removal Based on Space-Phase Method
Due to the large volume of geomagnetic data collected per unit area, it is not practical to employ standard outlier detection and removal/replacement methodologies.In a 20 m × 20 m survey area, an eight-sensor array configuration with 0.5 m separation may contain 150,000 data points collected at an average walking speed.This is an immense dataset in comparison to 1600 geomagnetic data points collected within a 20 m × 20 m survey grid with a 0.5 m transect and sampling spacing (Figure 5).Second, conventional outlier removal (despiking) methodologies mainly rely on a moving window where the operator expects a gridded data distribution.On the other hand, multi-sensor geomagnetic data usually has a complex spatial pattern.Therefore, in order to despike geomagnetic data, point patterns can be either regularized or different data processing algorithms can be used.In the IGEAN, outlier removal based on a space-phase method [14] was deployed after modifying the work by Mori [15].The phase-space of geomagnetic readings is based on u, the vertical magnetic gradient vector, where: Following Goring and Nikora [16], the postulate is that the magnitude of the geomagnetic gradient after a sensor visit along the transect should be less than a universal threshold.For an assumed normal random variable with a standard deviation of s, the expected absolute maximum of the threshold is calculated as: where is the universal threshold, n is the number of data points, and is the estimate of the standard deviation.The phase-space is enclosed by an ellipsoid with spherical coordinates p, f, and u [17]: Second, conventional outlier removal (despiking) methodologies mainly rely on a moving window where the operator expects a gridded data distribution.On the other hand, multi-sensor geomagnetic data usually has a complex spatial pattern.Therefore, in order to despike geomagnetic data, point patterns can be either regularized or different data processing algorithms can be used.In the IGEAN, outlier removal based on a space-phase method [14] was deployed after modifying the work by Mori [15].The phase-space of geomagnetic readings is based on u, the vertical magnetic gradient vector, where: Following Goring and Nikora [16], the postulate is that the magnitude of the geomagnetic gradient after a sensor visit along the transect should be less than a universal threshold.For an assumed normal random variable with a standard deviation of s, the expected absolute maximum of the threshold is calculated as: where γ u is the universal threshold, n is the number of data points, and σ is the estimate of the standard deviation.The phase-space is enclosed by an ellipsoid with spherical coordinates p, f, and u [17]: where a and b are major and minor axes of the ellipse projected on the u-∆ 2 u plane; c is the major axis in the ∆u-∆ 2 u; and a is the angle of rotation of the ellipsoid in the u-∆ 2 u plane.Geomagnetic readings falling outside of this ellipsoid are considered as data spikes (Figure 6).
Remote Sens. 2016, 8, 966 7 of 14 where a and b are major and minor axes of the ellipse projected on the u-Δ 2 u plane; c is the major axis in the Δu-Δ 2 u; and a is the angle of rotation of the ellipsoid in the u-Δ 2 u plane.Geomagnetic readings falling outside of this ellipsoid are considered as data spikes (Figure 6).

Data Sparsing based on Alpha Shapes
A multi-sensor array has its vast advantages, but resulting datasets also challenge the assumptions behind grid data-processing methodologies.In particular, many of the interpolation techniques perform better on regular data.Furthermore, off-the-shelf geophysical software mainly relies on the visualization of gridded data.Two main point-pattern issues have been observed in multi-sensor arrays that create concerns in the straightforward production of geomagnetic data maps (Figure 7a).
First, ground conditions or poor navigation may force data gaps between sensor arrays within the dense point pattern.Only good data collection practices can resolve this issue in the field and adjusting interpolation parameters to fill gaps should be avoided in most cases.Second, sensor arrays may overlap during navigation, resulting in denser data collection for those areas.Consequently, magnetic readings may differ for the same location as different sensor/array passes do not guarantee a fixed distance between the sensor and the ground or the same rotation angle between the sensor and the normal of the ground.In other words, the same location might bear two (or more) different geomagnetic readings in two (or more) successive passes of the sensor array from the same area.The best-practice followed during the IGEAN project was to keep the first collected geomagnetic reading at a location and remove the readings from other passes.
In order to establish this, a convex-hull algorithm roughly approximates the location and the extent of each transect.In return, a Boolean matrix is populated for documenting neighborhood topologies.The next step is the precise determination of overlapping areas between neighboring transects using alpha-shapes of each transect [18].This extra step is necessary since the convex-hull algorithm produces acceptable results only on convex shapes-a condition which is rarely satisfied for multi-array sensor transects.Once overlapping areas are established, secondary sensor readings are removed from the dataset so that an even distribution of dense point-patterns is satisfied (Figure 7b).

Data Sparsing Based on Alpha Shapes
A multi-sensor array has its vast advantages, but resulting datasets also challenge the assumptions behind grid data-processing methodologies.In particular, many of the interpolation techniques perform better on regular data.Furthermore, off-the-shelf geophysical software mainly relies on the visualization of gridded data.Two main point-pattern issues have been observed in multi-sensor arrays that create concerns in the straightforward production of geomagnetic data maps (Figure 7a).
First, ground conditions or poor navigation may force data gaps between sensor arrays within the dense point pattern.Only good data collection practices can resolve this issue in the field and adjusting interpolation parameters to fill gaps should be avoided in most cases.Second, sensor arrays may overlap during navigation, resulting in denser data collection for those areas.Consequently, magnetic readings may differ for the same location as different sensor/array passes do not guarantee a fixed distance between the sensor and the ground or the same rotation angle between the sensor and the normal of the ground.In other words, the same location might bear two (or more) different geomagnetic readings in two (or more) successive passes of the sensor array from the same area.The best-practice followed during the IGEAN project was to keep the first collected geomagnetic reading at a location and remove the readings from other passes.
In order to establish this, a convex-hull algorithm roughly approximates the location and the extent of each transect.In return, a Boolean matrix is populated for documenting neighborhood topologies.The next step is the precise determination of overlapping areas between neighboring transects using alpha-shapes of each transect [18].This extra step is necessary since the convex-hull algorithm produces acceptable results only on convex shapes-a condition which is rarely satisfied for multi-array sensor transects.Once overlapping areas are established, secondary sensor readings are removed from the dataset so that an even distribution of dense point-patterns is satisfied (Figure 7b).The alpha-shape approach can remove most of the data overlaps and even data distribution can be achieved for further processing.

Almyriotiki
Geomagnetic data from Magoula Almyriotiki [19] can be divided into three main interest groups (Figure 8).The first group is located at the top of the magoula, forming an arch with ~50 m radius.At least 34 architectural features fall into this category.Even though their shapes are not fully discernible it seems plausible to suggest that these buildings were mainly rectangular with areas ranging from 25-60 m 2 .The clear contrast between the geomagnetic values of "high positive-inside" and "low negative-outside" suggests pyro-related activities.Considering the spatial pattern of these anomalies, it is possible to speculate intentional firing.As another significant observation, there are a handful of anomalies at the outskirts of the magoula with the same geomagnetic signature of this group.If these structures were occupied simultaneously, then the dataset presents a curious case for the Neolithic intra-site settlement patterning.
The second group of anomalies are located at the lower skirts of the magoula.There are at least 55 structures in this group in which their rectangular shapes are immediately discernible.In some cases, internal divisions are visible, providing clues on the use of space.Considering negative geomagnetic measurements belonging to the "walls" of these features it can be suggested that diamagnetic materials, such as sedimentary stones were used in the construction of these architectural features.In another settlement, excavations at Koutroulou Magoula have revealed that limestone was used as the main material of building blocks of Neolithic houses [20].

Almyriotiki
Geomagnetic data from Magoula Almyriotiki [19] can be divided into three main interest groups (Figure 8).The first group is located at the top of the magoula, forming an arch with ~50 m radius.At least 34 architectural features fall into this category.Even though their shapes are not fully discernible it seems plausible to suggest that these buildings were mainly rectangular with areas ranging from 25-60 m 2 .The clear contrast between the geomagnetic values of "high positive-inside" and "low negative-outside" suggests pyro-related activities.Considering the spatial pattern of these anomalies, it is possible to speculate intentional firing.As another significant observation, there are a handful of anomalies at the outskirts of the magoula with the same geomagnetic signature of this group.If these structures were occupied simultaneously, then the dataset presents a curious case for the Neolithic intra-site settlement patterning.
The second group of anomalies are located at the lower skirts of the magoula.There are at least 55 structures in this group in which their rectangular shapes are immediately discernible.In some cases, internal divisions are visible, providing clues on the use of space.Considering negative geomagnetic measurements belonging to the "walls" of these features it can be suggested that diamagnetic materials, such as sedimentary stones were used in the construction of these architectural features.In another settlement, excavations at Koutroulou Magoula have revealed that limestone was used as the main material of building blocks of Neolithic houses [20].Geomagnetic prospection reveals further details about this group of features.Geomagnetic readings suggest a difference in post-depositional processes.Architectural features to the north present weaker magnetic signatures than the rest of the features in the same building assemblage.Assuming contemporaneity, it may be suggested that anomalies with weaker signatures are buried at larger depths within the ground or that the contrast between soils and structures is lower.This evidence comes in support of the impact of hydrological activity.The impact can be further observed on the settlement enclosure where a clearly defined anomaly at the eastern edges of the settlement is much harder to delineate to the north.
The third, and final, group of geomagnetism is related to a series of long linear anomalies which is more likely to be identified to settlement enclosures.Possibly, a ditch is documented to the east of the magoula, while the coverage to the west is not complete a similar anomaly is also observed, making a sharp 90 degrees turn to the east after running in north-south direction.Runnels et al. [21] report on two similar concentric ditches at Late Neolithic Makriyalos and suggest these features were the manifestations of Neolithic warfare.However, the amorphous shape of the enclosure to the north can be attributed to (paleo) hydrology around the site.This observation raises the possibility of ditches as preventive measures for flooding [22].

Perdika 1
Geomagnetic data from Perdika 1 [23] reveal a divergent image with respect to architectural elements at the site (Figure 9).The orientation of the buildings, and even some rooms, are visible.Though not as clear as the architecture of the settlement, there is evidence for a series of enclosures to the north of the site.Nevertheless, the immediate separation between the site and non-site provide information regarding the extent of the settlement.Geomagnetic prospection reveals further details about this group of features.Geomagnetic readings suggest a difference in post-depositional processes.Architectural features to the north present weaker magnetic signatures than the rest of the features in the same building assemblage.Assuming contemporaneity, it may be suggested that anomalies with weaker signatures are buried at larger depths within the ground or that the contrast between soils and structures is lower.This evidence comes in support of the impact of hydrological activity.The impact can be further observed on the settlement enclosure where a clearly defined anomaly at the eastern edges of the settlement is much harder to delineate to the north.
The third, and final, group of geomagnetism is related to a series of long linear anomalies which is more likely to be identified to settlement enclosures.Possibly, a ditch is documented to the east of the magoula, while the coverage to the west is not complete a similar anomaly is also observed, making a sharp 90 degrees turn to the east after running in north-south direction.Runnels et al. [21] report on two similar concentric ditches at Late Neolithic Makriyalos and suggest these features were the manifestations of Neolithic warfare.However, the amorphous shape of the enclosure to the north can be attributed to (paleo) hydrology around the site.This observation raises the possibility of ditches as preventive measures for flooding [22].

Perdika 1
Geomagnetic data from Perdika 1 [23] reveal a divergent image with respect to architectural elements at the site (Figure 9).The orientation of the buildings, and even some rooms, are visible.Though not as clear as the architecture of the settlement, there is evidence for a series of enclosures to the north of the site.Nevertheless, the immediate separation between the site and non-site provide information regarding the extent of the settlement.At Perdika 1 there is a clear difference of architectural practices at the site, depicted in two sets.The first type indicates the use of non-magnetic material as foundation material.Considering how the angles of rectangular buildings are preserved all over the site it is likely that the foundations were built out of slabs (e.g., limestone) or conglomerated pieces.This might, indeed, be the case as the settlement stands in close proximity to a creek as a potential source for non-magnetic building material.In this group of anomalies, architectural partitions are also visible in some cases.The second group is composed of buildings with high magnetic readings.It is possible that these buildings were made of mud-brick and burning resulted in distinct geomagnetic signatures.Furthermore, a cluster of buildings at the top of the magoula is indicated with high magnetic properties, suggesting a firing event.These buildings are tightly packed within a small space and differ from the rest of the features of the settlement in terms of their organizational layout.
A further area of interest is located to the north and to the east of the tightly packed cluster of buildings.The area is homogenous in terms of magnetic readings, suggesting an avoidance of use or frequent cleaning throughout occupation.Explaining the function and meaning of this space remains require further investigations, but remains important to show the diversity in the Neolithic use of space.Areas void of magnetic clusters are either located between architectural features and enclosures or between the elements of enclosures (see below).
The third feature of interest relates to the site enclosures.Even though the layout of these features is not complete, the existing evidence, especially to the north, is clear in terms of its morphology.It is likely that the settlement enclosure was fully encircling the settlement and surface taphonomic processes destroyed this feature in the expanse of protecting the inner settlement pattern.

Velestino-Mati
Velestino-Mati [24] is a bean/ellipsoid shaped settlement with two apparent settlement cores, one to the east and one to the west (Figure 10).Overall, the settlement is characterized by a series of buildings to the north and a sharp division between anthropogenic-natural divisions to the south.The core at the west contains more documented features than the eastern core despite the fact that some data are missing in the east.The saddle contains two anomalies circular in shape and high in magnetic values.Especially, the anomaly to the north might have been used as a large pit, serving both east and west settlement cores.At Perdika 1 there is a clear difference of architectural practices at the site, depicted in two sets.The first type indicates the use of non-magnetic material as foundation material.Considering how the angles of rectangular buildings are preserved all over the site it is likely that the foundations were built out of slabs (e.g., limestone) or conglomerated pieces.This might, indeed, be the case as the settlement stands in close proximity to a creek as a potential source for non-magnetic building material.In this group of anomalies, architectural partitions are also visible in some cases.The second group is composed of buildings with high magnetic readings.It is possible that these buildings were made of mud-brick and burning resulted in distinct geomagnetic signatures.Furthermore, a cluster of buildings at the top of the magoula is indicated with high magnetic properties, suggesting a firing event.These buildings are tightly packed within a small space and differ from the rest of the features of the settlement in terms of their organizational layout.
A further area of interest is located to the north and to the east of the tightly packed cluster of buildings.The area is homogenous in terms of magnetic readings, suggesting an avoidance of use or frequent cleaning throughout occupation.Explaining the function and meaning of this space remains require further investigations, but remains important to show the diversity in the Neolithic use of space.Areas void of magnetic clusters are either located between architectural features and enclosures or between the elements of enclosures (see below).
The third feature of interest relates to the site enclosures.Even though the layout of these features is not complete, the existing evidence, especially to the north, is clear in terms of its morphology.It is likely that the settlement enclosure was fully encircling the settlement and surface taphonomic processes destroyed this feature in the expanse of protecting the inner settlement pattern.

Velestino-Mati
Velestino-Mati [24] is a bean/ellipsoid shaped settlement with two apparent settlement cores, one to the east and one to the west (Figure 10).Overall, the settlement is characterized by a series of buildings to the north and a sharp division between anthropogenic-natural divisions to the south.The core at the west contains more documented features than the eastern core despite the fact that some data are missing in the east.The saddle contains two anomalies circular in shape and high in magnetic values.Especially, the anomaly to the north might have been used as a large pit, serving both east and west settlement cores.
The western core has features aligned in a circular shape, leaving an open space at the center.The architectural elements present evidence of burning.To the north of this circular layout, we also observe a series of potential buildings.It is possible that this alignment continues until the eastern core.However, chunks of data are missing to secure this argument.Even if one takes into account the missing data, the eastern core is not established as well as the western one.
Remote Sens. 2016, 8, 966 11 of 14 The western core has features aligned in a circular shape, leaving an open space at the center.The architectural elements present evidence of burning.To the north of this circular layout, we also observe a series of potential buildings.It is possible that this alignment continues until the eastern core.However, chunks of data are missing to secure this argument.Even if one takes into account the missing data, the eastern core is not established as well as the western one.
Figure 10.Geomagnetic data from Magoula Velestino-Mati.The bean shape settlement mound reveals some information about the built environment.There is little information on settlement enclosures but, when available, they are significantly distinct from the other two sites mentioned in this study.

Discussion
The geomagnetic prospecting under the IGEAN project explored Neolithic Thessaly in two interrelated domains: built environments and palaeolandscapes.Settlement enclosures with a variety of designs separate built environments from wider landscapes.The syntheses of results were possible only after acquiring the complete coverage of sites with geomagnetic prospection, following Kvamme [2] in his geophysical surveys as a landscape archaeology approach.

Built Environments
The geophysical data shows that there is a significant variety in structural organization.These results trigger a new frontier of questions.For instance, the results from Perdika 1 revealed a rectangular layout in comparison to the circular arrangement of mounded settlements.IGEAN suggests that almost all tell-sites differ in terms of their circularity.More questions arise regarding the existence of clusters/neighborhoods (Velestino-Mati) or the differences in house sizes and their internal organizations (Almyriotiki vs. Perdika 1) that may directly have been associated with the social organization in Neolithic society.Finally, the geophysical results indicate the first evidence of settlement expansion giving clues on the growth of communities.

Enclosures
Where enclosures have been detected, they are generally identified as single, double, or multiple ditch-like features which have a large variety of sizes and shapes.Taking into account the lack of knowledge about the function of these enclosures, it was possible to investigate the diversity of these features in two ways: (1) by examining if specific enclosure forms are related to the type of settlement layout; and (2) by investigating if there is an association between enclosures and their immediate fluvial system.Considering the hypothetic practice of floodwater farming strategy in Neolithic Figure 10.Geomagnetic data from Magoula Velestino-Mati.The bean shape settlement mound reveals some information about the built environment.There is little information on settlement enclosures but, when available, they are significantly distinct from the other two sites mentioned in this study.

Discussion
The geomagnetic prospecting under the IGEAN project explored Neolithic Thessaly in two inter-related domains: built environments and palaeolandscapes.Settlement enclosures with a variety of designs separate built environments from wider landscapes.The syntheses of results were possible only after acquiring the complete coverage of sites with geomagnetic prospection, following Kvamme [2] in his geophysical surveys as a landscape archaeology approach.

Built Environments
The geophysical data shows that there is a significant variety in structural organization.These results trigger a new frontier of questions.For instance, the results from Perdika 1 revealed a rectangular layout in comparison to the circular arrangement of mounded settlements.IGEAN suggests that almost all tell-sites differ in terms of their circularity.More questions arise regarding the existence of clusters/neighborhoods (Velestino-Mati) or the differences in house sizes and their internal organizations (Almyriotiki vs. Perdika 1) that may directly have been associated with the social organization in Neolithic society.Finally, the geophysical results indicate the first evidence of settlement expansion giving clues on the growth of communities.

Enclosures
Where enclosures have been detected, they are generally identified as single, double, or multiple ditch-like features which have a large variety of sizes and shapes.Taking into account the lack of knowledge about the function of these enclosures, it was possible to investigate the diversity of these features in two ways: (1) by examining if specific enclosure forms are related to the type of settlement layout; and (2) by investigating if there is an association between enclosures and their immediate fluvial system.Considering the hypothetic practice of floodwater farming strategy in Neolithic Thessaly [22], IGEAN made further investigations to understand if enclosures are actually related to paleohydrological activity around the settlements.

Paleoenvironments
The IGEAN project also investigated areas surrounding Neolithic settlements.First, the IGEAN project was able to locate numerous palaeochannels around Neolithic settlements.Considering the importance of hydrological systems in the livelihood of early farming communities it becomes clear that examining whether flooding and related stream flow phenomena is critical to understand.A recorded rate of subsidence of 1.5 m/1000 years for the Larissa basin has been reported, signifying the necessity of investigating the flooding susceptibility of this particular region [25].
In some cases, (e.g., Almyriotiki) it was clear that past flooding activities adversely affected the preservation of settlement ditches.Even though it is not possible to date the exact dating of this flooding event, the phenomenon immediately raises the relationship between flooding and the function of these ditches.This comes as an alternative hypothesis for enclosures as defensive features as suggested by Runnels et al. [21].The same phenomenon can be also argued for the enclosures of Perdika 1, albeit with less secure evidence of a flooding event.

Conclusions
The IGEAN project served as a prime example for revealing the efficiency of multi-sensor geomagnetic prospecting.The complete data coverage of Neolithic settlements enables researchers to pose key questions; this has not been possible before the advent of array systems in archaeological research.Exceptions are only in the expanse of enormous time spent on data collection and analyses [26].The power of multi-sensor geomagnetic prospecting is evident in its wide coverage in a short period of time.Array magnetometry is invaluable for revealing most of the significant features in the landscape, but accuracy of interpretations drastically increase when other methods are also deployed.In such an effort Simon et al. [27] discuss the geometry of enclosures at Almyriotiki using the results of ground penetrating radar (GPR) and electromagnetic induction (EM) surveys (Figure 11).Similarly, the accurate characterization of a large building encountered at the same site is only possible when GPR is included in the survey schema (Figure 12).Despite these illustrations, until the time comes for the efficient multi-sensor instrumentation of other geophysical techniques, (multi-sensor) magnetometry-once again-remains as the workhorse of archaeological prospecting [28].
Remote Sens. 2016, 8, 966 12 of 14 Thessaly [22], IGEAN made further investigations to understand if enclosures are actually related to paleohydrological activity around the settlements.

Paleoenvironments
The IGEAN project also investigated areas surrounding Neolithic settlements.First, the IGEAN project was able to locate numerous palaeochannels around Neolithic settlements.Considering the importance of hydrological systems in the livelihood of early farming communities it becomes clear that examining whether flooding and related stream flow phenomena is critical to understand.A recorded rate of subsidence of 1.5 m/1000 years for the Larissa basin has been reported, signifying the necessity of investigating the flooding susceptibility of this particular region [25].
In some cases, (e.g., Almyriotiki) it was clear that past flooding activities adversely affected the preservation of settlement ditches.Even though it is not possible to date the exact dating of this flooding event, the phenomenon immediately raises the relationship between flooding and the function of these ditches.This comes as an alternative hypothesis for enclosures as defensive features as suggested by Runnels et al [21].The same phenomenon can be also argued for the enclosures of Perdika 1, albeit with less secure evidence of a flooding event.

Conclusions
The IGEAN project served as a prime example for revealing the efficiency of multi-sensor geomagnetic prospecting.The complete data coverage of Neolithic settlements enables researchers to pose key questions; this has not been possible before the advent of array systems in archaeological research.Exceptions are only in the expanse of enormous time spent on data collection and analyses [26].The power of multi-sensor geomagnetic prospecting is evident in its wide coverage in a short period of time.Array magnetometry is invaluable for revealing most of the significant features in the landscape, but accuracy of interpretations drastically increase when other methods are also deployed.In such an effort Simon et al. [27] discuss the geometry of enclosures at Almyriotiki using the results of ground penetrating radar (GPR) and electromagnetic induction (EM) surveys (Figure 11).Similarly, the accurate characterization of a large building encountered at the same site is only possible when GPR is included in the survey schema (Figure 12).Despite these illustrations, until the time comes for the efficient multi-sensor instrumentation of other geophysical techniques, (multisensor) magnetometry-once again-remains as the workhorse of archaeological prospecting [28].

Figure 1 .
Figure 1.Thessaly was the gateway for the Neolithization of Europe.The Innovative Geophysical Approaches for the Study of Early Agricultural Villages of Neolithic Thessaly (IGEAN) project examined 24 settlements in detail using geospatial technologies.Three sites (Almyriotiki, Perdika 1, and Velestino-Mati) constitute the demonstration case studies of this paper.Within a new scientific agenda of investigating Neolithic cultures, the Laboratory of Geophysical-Satellite Remote Sensing and Archaeo-environment (GeoSat ReSeArch Lab) of the Foundation for Research and Technology-Hellas (FORTH) initiated (2013-2015) a remote sensing campaign to study the physical landscape and dynamics of Neolithic settlements within the coastal hinterlands of Eastern Thessaly.The IGEAN (http://igean.ims.forth.gr/)research project is the first systematic and extensive geophysical survey of Neolithic magoules in the region and is the first to utilize multi-sensor geomagnetic sensors in the Eastern Mediterranean.This paper focuses on the geomagnetic prospection component of the IGEAN and reveals three prime examples from the project.In doing so, it aims to show the power of multi-sensor systems as an emerging paradigm in archaeology.Given that the surface conditions are appropriate for extensive data collection, the efficiency of multi-sensors over single or dual sensor systems is beyond doubt.It is now truly possible to conduct "geophysical surveys as landscape archaeology"[2].

Figure 1 .
Figure 1.Thessaly was the gateway for the Neolithization of Europe.The Innovative Geophysical Approaches for the Study of Early Agricultural Villages of Neolithic Thessaly (IGEAN) project examined 24 settlements in detail using geospatial technologies.Three sites (Almyriotiki, Perdika 1, and Velestino-Mati) constitute the demonstration case studies of this paper.

Figure 2 .
Figure 2. The SENSYS Sensorik and Systemtechnologie can carry up to 16 sensors, offering variations for sensor numbers and separations.The IGEAN setup included eight sensors with 50 cm separation.A high-precision Global Navigation Satellite System (GNSS) receiver helps in navigation and sensor-data registration.

Figure 3 .
Figure 3.The IGEAN workflow for processing geomagnetic data.The first three steps (ASCII parsing, quality control, and data clipping) are standardized.The proceeding steps and their related parameters are case dependent.

Figure 4 .
Figure 4.The raw data (above) is imported from the SENSYS software in GeoTIFF format and displayed in a Geographic Information System (GIS).The processed data (below) shows drastic improvements over the original dataset.The data is an excerpt from Magoula Almyriotiki (see Section 5.1).

Figure 3 .
Figure 3.The IGEAN workflow for processing geomagnetic data.The first three steps (ASCII parsing, quality control, and data clipping) are standardized.The proceeding steps and their related parameters are case dependent.

Figure 3 .
Figure 3.The IGEAN workflow for processing geomagnetic data.The first three steps (ASCII parsing, quality control, and data clipping) are standardized.The proceeding steps and their related parameters are case dependent.

Figure 4 .
Figure 4.The raw data (above) is imported from the SENSYS software in GeoTIFF format and displayed in a Geographic Information System (GIS).The processed data (below) shows drastic improvements over the original dataset.The data is an excerpt from Magoula Almyriotiki (see Section 5.1).

Figure 4 .
Figure 4.The raw data (above) is imported from the SENSYS software in GeoTIFF format and displayed in a Geographic Information System (GIS).The processed data (below) shows drastic improvements over the original dataset.The data is an excerpt from Magoula Almyriotiki (see Section 5.1).

Figure 5 .
Figure 5.The grid collection of geomagnetic data (red triangles) results in far fewer data points than multi-sensor collection.The comparison is the most clear in the figure inlay.For both images, the north is up.

Figure 5 .
Figure 5.The grid collection of geomagnetic data (red triangles) results in far fewer data points than multi-sensor collection.The comparison is the most clear in the figure inlay.For both images, the north is up.

Figure 6 .
Figure 6.The space-phase method can detect data outliers in an efficient manner.Geomagnetic readings falling outside of the 3D ellipsoid (in red) are tagged as data spikes.

Figure 6 .
Figure 6.The space-phase method can detect data outliers in an efficient manner.Geomagnetic readings falling outside of the 3D ellipsoid (in red) are tagged as data spikes.

Figure 7 .
Figure 7. Distribution of data points is determined by track navigation.(a) Horizontal hatching indicates overlapping areas and thus these areas contain denser data.Diagonal hatching indicates data gaps, which can be avoided by good navigation skills in most cases.(b) The alpha-shape approach can remove most of the data overlaps and even data distribution can be achieved for further processing.

Figure 7 .
Figure 7. Distribution of data points is determined by track navigation.(a) Horizontal hatching indicates overlapping areas and thus these areas contain denser data.Diagonal hatching indicates data gaps, which can be avoided by good navigation skills in most cases; (b) The alpha-shape approach can remove most of the data overlaps and even data distribution can be achieved for further processing.

Figure 8 .
Figure 8. Geomagnetic data from Magoula Almyriotiki.Two distinct architectural groups and settlement enclosures are clearly visible.

Figure 8 .
Figure 8. Geomagnetic data from Magoula Almyriotiki.Two distinct architectural groups and settlement enclosures are clearly visible.

Figure 9 .
Figure 9. Geomagnetic data from Perdika 1. Palimpsest of two architectural groups are visible in the dataset.The settlement enclosure is the most visible in the north.

Figure 9 .
Figure 9. Geomagnetic data from Perdika 1. Palimpsest of two architectural groups are visible in the dataset.The settlement enclosure is the most visible in the north.

Figure 11 .
Figure 11.Magnetic susceptibility data from the enclosure in Magoula Almyriotiki suggest a single running feature; unlike the geomagnetic data which reveal a double-ditch like feature.

Figure 11 .
Figure 11.Magnetic susceptibility data from the enclosure in Magoula Almyriotiki suggest a single running feature; unlike the geomagnetic data which reveal a double-ditch like feature.

Figure 12 .
Figure 12.Magnetometers are the work-horses of archaeological prospection.However, it is the best practice to integrate other sensors into the work schema.(a) Geomagnetic data from Magoula Almyriotiki reveals a large building to the south of the settlement.(b) GPR survey in the same area further corrects this argument and suggests three buildings; two rectangular buildings to the west and another one, equally partitioned, to the east.

Figure 12 .
Figure 12.Magnetometers are the work-horses of archaeological prospection.However, it is the best practice to integrate other sensors into the work schema.(a) Geomagnetic data from Magoula Almyriotiki reveals a large building to the south of the settlement; (b) GPR survey in the same area further corrects this argument and suggests three buildings; two rectangular buildings to the west and another one, equally partitioned, to the east.

Remote Sens. 2016, 8, 966 4 of 14 Figure 2. The
SENSYS Sensorik and Systemtechnologie can carry up to 16 sensors, offering variations for sensor numbers and separations.The IGEAN setup included eight sensors with 50 cm separation.A high-precision Global Navigation Satellite System (GNSS) receiver helps in navigation and sensordata registration.