Integration of Remote-Sensing Techniques for the Preventive Conservation of Paleolithic Cave Art in the Karst of the Altamira Cave

: Rock art offers traces of our most remote past and was made with mineral and organic substances in shelters, walls, or the ceilings of caves. As it is notably fragile, it is fortunate that some instances remain intact—but a variety of natural and anthropogenic factors can lead to its disappearance. Therefore, as a valuable cultural heritage, rock art requires special conservation and protection measures. Geomatic remote-sensing technologies such as 3D terrestrial laser scanning (3DTLS), drone ﬂight, and ground-penetrating radar (GPR) allow us to generate exhaustive documentation of caves and their environment in 2D, 2.5D, and 3D. However, only its combined use with 3D geographic information systems (GIS) lets us generate new cave maps with details such as overlying layer thickness, sinkholes, fractures, joints, and detachments that also more precisely reveal interior–exterior interconnections and gaseous exchange; i.e., the state of senescence of the karst that houses the cave. Information of this kind is of great value for the research, management, conservation, monitoring, and dissemination of cave art.

From a remote-sensing viewpoint, caves are characterized by particular conditions, such as complex access, confined spaces, high humidity, constant temperature, absence of natural light, a potential relationship between the interior and the exterior, etc. [4,21].
Separately, photogrammetric unmanned aerial vehicle (UAV) flights can provide accurate terrain models [22,23], 3D terrestrial laser scanner (3DTLS) surveys provide detailed and correct cave models, while ground-penetrating radar (GPR) provides information on the particularities of the rock mass under investigation.
Separately, photogrammetric unmanned aerial vehicle (UAV) flights can provide accurate terrain models [22,23], 3D terrestrial laser scanner (3DTLS) surveys provide detailed and correct cave models, while ground-penetrating radar (GPR) provides information on the particularities of the rock mass under investigation.
It is the combined use of these tools that enables the generation of new cartographic information on the features of cave karst complexes, such as rock thickness, fractures, joints, sinkholes, cave seepage, water dripping, stability analysis, and gaseous exchanges, which are valuable for conservation studies of the cave.
The Altamira Cave is an exceptional cultural heritage site that is internationally recognized as a masterpiece of History of Art, being added to the UNESCO World Heritage List in 1985.It is in Santillana del Mar (Cantabria) and is 296 m long, widening in some areas, such as the Polychrome Hall.It is the first cave in which the existence of rock art from the Palaeolithic period was identified (Figure 1).The cave was discovered in 1868, while the cave art was discovered in 1879 by Marcelino Sanz de Sautuola and his daughter María.The authenticity of the cave paintings The cave was discovered in 1868, while the cave art was discovered in 1879 by Marcelino Sanz de Sautuola and his daughter María.The authenticity of the cave paintings was accepted by the scientific community in 1902 [24,25].It has paintings and engravings that were made from over 35,000 years ago to 14,000 years ago, among which the set of polychrome bison stands out.
Since its discovery, it has been the object of special conservational attention and care.It is in the upper part of an ancient karst (Figure 2a) that has a tendency to disappear by collapse and subsidence: due to the tabular structure of the limestone and the steep slope of the fracture planes, the geological evolution of Altamira is more marked by a series of gravitational collapses of the limestone rock layers than by the circulation of water [26].
ronmental conservation problems that continue to this day.Subsequent interventions focused on trying to reverse the infiltration of rainwater through the network of fissures and fractures in the cave, which is the main problem for the conservation of the paintings.Artificial walls were built inside the cave (Figure 2a), some cracks were filled with hydraulic mortar, and a layer of Portland cement was spread over the exterior surface [27].
The current state of conservation of Altamira results from structural geological problems due to the accumulation over time of multiple transformations and micro-environmental disturbances.Refurbishment works, changes in the external and internal soils, archaeological excavations, mass visits, and other activities have generated conservation conditions different from those that the cave had before its discovery, where its substantial conservation was caused by its low water infiltration rate and its stable microclimatic conditions, the latter being caused by its natural closure [28].The first writings in the local and national press reported the existence of collapses throughout the cave.Faced with this structural instability in some parts of the cave at the beginning of the 20th century, important transformations were proposed, such as the raising, in the vestibular area, of a large stone wall to act as a support for the strata of the ceiling, which was threatening to collapse.This intervention, which definitively isolated the Polychrome room from the rest of the cavity, has generated some of the microenvironmental conservation problems that continue to this day.Subsequent interventions focused on trying to reverse the infiltration of rainwater through the network of fissures and fractures in the cave, which is the main problem for the conservation of the paintings.Artificial walls were built inside the cave (Figure 2a), some cracks were filled with hydraulic mortar, and a layer of Portland cement was spread over the exterior surface [27].
The current state of conservation of Altamira results from structural geological problems due to the accumulation over time of multiple transformations and micro-environmental disturbances.Refurbishment works, changes in the external and internal soils, archaeological excavations, mass visits, and other activities have generated conservation conditions different from those that the cave had before its discovery, where its substantial conservation was caused by its low water infiltration rate and its stable microclimatic conditions, the latter being caused by its natural closure [28].
The most worrying factors in the conservation of the Altamira cave concern deterioration processes related to geological structural stability and infiltration water, as well as biodeterioration and human presence.To understand their dynamics and to minimize the alterations generated by them, it is necessary to develop effective conservation strategies aimed at controlling the different agents of deterioration.The water that enters the cave through the network of fissures and fractures in the overlying layer plays an important role in the processes of dissolution and dragging of the supporting rock.The processes of microcorrosion, loss of adhesion, cracking, washing, and loss of support are worrying and are still active today.
Knowing and studying the underlying dynamics of these processes hypothetically associated with the destruction of paint at Altamira is essential to design strategies to control each alteration agent and the interactions between them.
The aim of this research was to improve the understanding of the karst system surrounding Altamira Cave by combining remote-sensing technologies for its preventive conservation, thus contributing to decision-making for its management.
An exhaustive and combined study of remote-sensing techniques was carried out, after which we performed 3D cartography to integrate the exterior, the inner topography of the cave, and the main discontinuities and karst features detected around.As we will show, this 3D integrated mapping is a crucial tool for monitoring the main fluid flows into the cave and other conservation purposes.
In the paper, Section 2 shows the data sources integrated to carry out the study, i.e., GNSS in different operational modes, 3DTLS for the geometrical characterization of the cave, a UAV flight to model the terrain, and GPR for the overlying layer.Each subsection describes the accuracies achieved with each method.Section 3 shows the results obtained from integrating two or more techniques and how to visualize them.Finally, Section 4 contains the discussion.

Materials and Methods
Caves with rock art are formed by the dissolution of limestone rock by slightly acidic water, which may involve processes of chemical and geological activity, tectonic forces, and atmospheric influences.They are irregular in shape and variable in extent, so it is necessary to integrate different georeferenced remote-sensing techniques to derive new cartographies and generate new information on the cultural heritage element [29].
The flowchart that follows first required the creation of a reference frame, which was observed using the static mode of GNSS, obtaining an average error of 1.7 cm at the vertices.This reference frame was observed and adjusted as a microgeodetic network, obtaining an average accuracy of 1.1 cm in the determination of vertex coordinates.
The UAV flight was georeferenced by GNSS in RTK mode.Fifty-two control points were taken throughout the area and 40 points were taken to validate the model.The result of the validation points was an error of 2.1 cm in planimetry and 2.8 cm in elevation.
The GPR grids were staked using GNSS in RTK mode.To keep the direction constant, a canvas was attached to each grid.From these efforts, the three models were generated and then integrated to be analyzed as a whole (Figure 3).

Global Navigation Satellite System (GNSS)
GNSS (Global Navigation Satellite System) refers to the set of satellite navigation system technologies that provide geospatial positioning with autonomous global coverage.
The origins of GNSS date to the 1970s with the development of the US military GPS (Global Positioning System), which, through a network of satellites, a GNSS receiver can determine its position in four dimensions (longitude, latitude, altitude, and time) [30].
In September 2013, a network of points was created using three TOPCON Hyper II [31] receivers simultaneously, 2 of which were moved to other vertices after 1 h of observation to create the reference frame, which was the beginning of the traverse of the cave.The geodetic reference system used to create the reference frame was the European Terrestrial Reference System 1989 (ETRS89).
The software Topcon Tools [32] was used to calculate observed phase changes to form the differencing observation equations, perform the least-squares adjustment, and output the adjusted baseline vector components.The mean accuracy of the reference frame was 1.7 cm.
GPR grids were staked out by using the Real Time Kinematic (RTK) mode.The RTK positioning technique is based on the carrier solution of the signals transmitted by the global navigation satellite systems GPS; in our case, Glonass and Navstar.A reference station provides instantaneous corrections for mobile stations, which brings the obtained accuracy to the centimeter level.The base station retransmits the phase of the carrier it measured, and the mobile units compare their own phase measurements with that received from the reference station.This lets the mobile stations calculate their relative positions with millimeter accuracy, with their absolute relative positions being related to the base station coordinates.This technique requires the availability of at least one reference station with known coordinates that has a GNSS receiver and a radio transmitter modem.The station generates and transmits differential corrections for the stations, which use the data to precisely determine their positions.[33,34]

3D Terrestrial Laser Scanning
Having an accurate cartographic model of a cave is important when knowledgebased prediction models need to be estimated from georeferenced parameters such as microbiological, climatic, hydrochemical, geomorphological, etc. that are highly correlated either between them or with what happens in the exterior.[18].
Starting from a vertex of the reference frame "b-2" (Figure 2b), a closed traverse backand-forth was made that was adjusted and compensated through the cave using a TOP-CON GPT-7503 total station [35].Sixteen traverse stations were needed, and the total

Global Navigation Satellite System (GNSS)
GNSS (Global Navigation Satellite System) refers to the set of satellite navigation system technologies that provide geospatial positioning with autonomous global coverage.
The origins of GNSS date to the 1970s with the development of the US military GPS (Global Positioning System), which, through a network of satellites, a GNSS receiver can determine its position in four dimensions (longitude, latitude, altitude, and time) [30].
In September 2013, a network of points was created using three TOPCON Hyper II [31] receivers simultaneously, 2 of which were moved to other vertices after 1 h of observation to create the reference frame, which was the beginning of the traverse of the cave.The geodetic reference system used to create the reference frame was the European Terrestrial Reference System 1989 (ETRS89).
The software Topcon Tools [32] was used to calculate observed phase changes to form the differencing observation equations, perform the least-squares adjustment, and output the adjusted baseline vector components.The mean accuracy of the reference frame was 1.7 cm.
GPR grids were staked out by using the Real Time Kinematic (RTK) mode.The RTK positioning technique is based on the carrier solution of the signals transmitted by the global navigation satellite systems GPS; in our case, Glonass and Navstar.A reference station provides instantaneous corrections for mobile stations, which brings the obtained accuracy to the centimeter level.The base station retransmits the phase of the carrier it measured, and the mobile units compare their own phase measurements with that received from the reference station.This lets the mobile stations calculate their relative positions with millimeter accuracy, with their absolute relative positions being related to the base station coordinates.This technique requires the availability of at least one reference station with known coordinates that has a GNSS receiver and a radio transmitter modem.The station generates and transmits differential corrections for the stations, which use the data to precisely determine their positions [33,34].

3D Terrestrial Laser Scanning
Having an accurate cartographic model of a cave is important when knowledgebased prediction models need to be estimated from georeferenced parameters such as microbiological, climatic, hydrochemical, geomorphological, etc. that are highly correlated either between them or with what happens in the exterior [18].
Starting from a vertex of the reference frame "b-2" (Figure 2b), a closed traverse back-and-forth was made that was adjusted and compensated through the cave using a TOPCON GPT-7503 total station [35].Sixteen traverse stations were needed, and the total length was 430 m.The angular closure error was 0.0218 g , and the linear closure error was X = −0.001m, Y = −0.005m, Z = 0.A total of 66 checkerboards distributed around the cave were radiated from the vertices of the traverse; these are the reference system of the 3D laser scanner, on which the scans are georeferenced.(Figure 2b).
The field campaign was carried out in December 2013; a FARO FOCUS X-130 [36] was used, and 300 scans were required.Calibrated spheres were used as tie points and the traverse checkerboards as references, which were adjusted with an accuracy of 2.7 mm for 95% of the points (Figure 2a) [37].
This new model updated the traditional values of the Altamira cave, which are as follows: • Length: 296.5 m (20 longer) The 3D model enabled us to obtain, apart from the classical cartographies of plant and section planes, other cartographies such as the height of roofs, floors, and galleries.

Unmanned Aerial Vehicle
Currently, the use of drones for surveying and mapping work is widespread.Drones allow surveys of large areas to be carried out with high reliability in a shorter time than using traditional techniques such as Topographic Total Stations (TTS) and GNSS.UAVs make it possible to obtain something fundamental for our model, which is to have continuous, sufficiently precise information of the entire work area, including places that are difficult to access: in our case, the sinkholes [38][39][40].
The platform selected for the execution of the photogrammetric flights was a TOPCON Intel Falcon 8+ Drone [41], on which a Sony A7 R Mark ii full-frame camera [42] was mounted with a Sony Sonnar T* FE 35 mm F/2.8 lens.The flights were carried out with the flight-planning system in a simple and fast way.
To reach a ground sample distance (GSD) of 2 cm, the flight-planning system was used to obtain measurements throughout the entire cave environment (40 Ha).More than 4800 pictures were taken, 52 Ground Control Points (GCP) were measured with RTK-GNSS to reference the flight, and another 40 points were measured for quality control.The mean error of the control points was 1.76 cm.
This allowed us to obtain 2D information such as orthoimages, 2.5D information such as a Digital Surface Model, which was filtered to remove trees to obtain a Digital Elevation Model (DEM), and 3D information such as point clouds and a 3D model of the whole area.The 3D model of the cave environment (40 Ha) was created at a definition of 12 million triangles textured with 8 textures of 8192 × 8192 pixels (Figure 4), while the model of the area next to the cave (4 Ha) was created with 4 million polygons and 4 textures of 8192 × 8192 pixels, which is the model we used for integrating the rest of the data.Both models were exported to .Obj, a geometry definition file format, to be integrated with other data.

Ground-Penetrating Radar
Ground-Penetrating Radar (GPR) is commonly used in the investigation of the internal structure of rock masses, which are generally defined as mechanically discontinuous, anisotropic and heterogeneous environments [43][44][45].GPR provides continuous, high-resolution information on the internal structure of a rock mass and its physical properties in both vertical and lateral directions [46][47][48][49].The data recorded with the GPR (radargrams) depend on the transmission of the electromagnetic wave in the medium and the amount of energy reflected, diffracted, or critically refracted by the medium.Therefore, the propagation of a wave depends on the electromagnetic properties (magnetic permittivity, electrical conductivity, and dielectric permittivity) of the medium through which it is transmitted and the contrast of these electromagnetic properties with those of other media.

Field Data Acquisition
It was necessary to collect data from the whole site of the Altamira Cave and its surroundings.A three-dimensional (3D) GPR methodology was carried out to pinpoint and define the existence and geometry of the main discontinuities and fractures in the karst system where the cave is located.One feature of GPR surveying is the inevitable decrease in resolution with depth.The frequency of the antenna is a crucial choice.Different resolutions and penetration depths can be achieved using antennas with different nominal center frequencies: higher resolution and shallower penetration depths are typically obtained with higher frequency antennas, while deeper penetrations and lower resolution are achieved with lower frequency antennas [61,62].A commercial GPR equipment SIR-3000 (GSSI) was used for the study.In this study, we present the results of three GPR grids.

Ground-Penetrating Radar
Ground-Penetrating Radar (GPR) is commonly used in the investigation of the internal structure of rock masses, which are generally defined as mechanically discontinuous, anisotropic and heterogeneous environments [43][44][45].GPR provides continuous, highresolution information on the internal structure of a rock mass and its physical properties in both vertical and lateral directions [46][47][48][49].The data recorded with the GPR (radargrams) depend on the transmission of the electromagnetic wave in the medium and the amount of energy reflected, diffracted, or critically refracted by the medium.Therefore, the propagation of a wave depends on the electromagnetic properties (magnetic permittivity, electrical conductivity, and dielectric permittivity) of the medium through which it is transmitted and the contrast of these electromagnetic properties with those of other media.

Field Data Acquisition
It was necessary to collect data from the whole site of the Altamira Cave and its surroundings.A three-dimensional (3D) GPR methodology was carried out to pinpoint and define the existence and geometry of the main discontinuities and fractures in the karst system where the cave is located.One feature of GPR surveying is the inevitable decrease in resolution with depth.The frequency of the antenna is a crucial choice.Different resolutions and penetration depths can be achieved using antennas with different nominal center frequencies: higher resolution and shallower penetration depths are typically obtained with higher frequency antennas, while deeper penetrations and lower resolution are achieved with lower frequency antennas [61,62].A commercial GPR equipment SIR-3000 (GSSI) was used for the study.In this study, we present the results of three GPR grids. - In order to pinpoint discontinuities and karst structures in deep layers, a 100 MHz nominal center frequency antenna was used for the area in and around Altamira Cave.The 100 MHz antenna was able to penetrate to depths of about 10 m in the study area, and a centimeter resolution could be achieved [63,64].The design of the GPR profile grid mainly took into account the dimensions and directionalities of a series of karst structures and faults, in some cases extending over more than half a kilometer, mapped in previous studies in the Altamira Cave area [65,66].Therefore, a profile grid for GPR data acquisition was designed with equidistant 2D orthogonal lines, spaced every 3 m, with respect to a longitudinal axis that was, in turn, perpendicular to the main fracture directions in the Altamira Cave area, as described in these previous studies.Every 2D profile line was topographically georeferenced by GNSS (Figures 5 and 6).
Remote Sens. 2023, 15, x FOR PEER REVIEW 8 of 25 -In order to pinpoint discontinuities and karst structures in deep layers, a 100 MHz nominal center frequency antenna was used for the area in and around Altamira Cave.The 100 MHz antenna was able to penetrate to depths of about 10 m in the study area, and a centimeter resolution could be achieved [63,64].The design of the GPR profile grid mainly took into account the dimensions and directionalities of a series of karst structures and faults, in some cases extending over more than half a kilometer, mapped in previous studies in the Altamira Cave area [65,66].Therefore, a profile grid for GPR data acquisition was designed with equidistant 2D orthogonal lines, spaced every 3 m, with respect to a longitudinal axis that was, in turn, perpendicular to the main fracture directions in the Altamira Cave area, as described in these previous studies.Every 2D profile line was topographically georeferenced by GNSS (Figures 5 and 6).-Due to the importance of the cave paintings on the Polychrome Hall ceiling for the Cave of Altamira, it was considered necessary to study the first meters of its overlying rock strata and its lapiés area in greater detail and at a higher resolution.Two GPR profile grids were planned in the Polychrome Hall area: o A profile grid for the study of its strata in the first meters using a 400 MHz nominal center frequency antenna to identify discontinuities and karst structures.The 400 MHz antenna was able to penetrate to depths of about 4 m in the study area, with centimeter resolution [67].o A profile grid for the study of its lapiés zone using a 900 MHz nominal center frequency antenna to define the lapiés groove system.The 900 MHz antenna was able to penetrate to depths of about 1.5 m in the study area with a centimeter resolution [68].
These two grids were designed according to the extensive system of lineaments in the exterior and interior of Altamira Cave, mainly in the form of joints and with predominant directions N45, N68, N110, N168, and N-S [69].Accordingly, both profile grids for GPR data acquisition were designed with equidistant 2D orthogonal lines, spaced every 1 m, with respect to a longitudinal axis that was, in turn, perpendicular to the main joints directions in the Polychrome Hall area, as described in the aforementioned previous study.Every 2D profile line was topographically georeferenced by GNSS (Figures 7 and  8).-Due to the importance of the cave paintings on the Polychrome Hall ceiling for the Cave of Altamira, it was considered necessary to study the first meters of its overlying rock strata and its lapiés area in greater detail and at a higher resolution.Two GPR profile grids were planned in the Polychrome Hall area: • A profile grid for the study of its strata in the first meters using a 400 MHz nominal center frequency antenna to identify discontinuities and karst structures.The 400 MHz antenna was able to penetrate to depths of about 4 m in the study area, with centimeter resolution [67].

•
A profile grid for the study of its lapiés zone using a 900 MHz nominal center frequency antenna to define the lapiés groove system.The 900 MHz antenna was able to penetrate to depths of about 1.5 m in the study area with a centimeter resolution [68].
These two grids were designed according to the extensive system of lineaments in the exterior and interior of Altamira Cave, mainly in the form of joints and with predominant directions N45, N68, N110, N168, and N-S [69].Accordingly, both profile grids for GPR data acquisition were designed with equidistant 2D orthogonal lines, spaced every 1 m, with respect to a longitudinal axis that was, in turn, perpendicular to the main joints directions in the Polychrome Hall area, as described in the aforementioned previous study.Every 2D profile line was topographically georeferenced by GNSS (Figures 7 and 8).In the field campaign, 454 profiles were recorded: 329 profiles totalling 20,850 m 2 were surveyed with the 100 MHz antenna, 76 profiles totalling 850 m 2 with the 400 MHz antenna, and 49 profiles totalling 550 m 2 with the 900 MHz antenna.

Data Processing Velocity Estimation
The GPR reflection data was acquired in the same week under very similar weather conditions with no rainfall.Thus, the subsurface moisture conditions in this investigated area could be considered almost constant for the investigated depth during the acquisition of the GPR field data.This velocity estimation assumption was made for the velocity-todepth conversions, considering that possible changes in subsurface conditions could cause velocity variations that would affect these velocity-to-depth conversions [61,69].
The average velocity of the GPR wave was determined by the hyperbola fitting method on a set of hyperbolas recorded in different profiles, obtaining an average velocity value of 11.17 cm/ns.The dielectric permittivity (ɛ) was calculated to be 7.2, according to the following equation for low-loss media [61,62]: where v is the electromagnetic wave velocity and c is the velocity of light in free space (c ≈ 29.97 cm/ns).This value was applied for the calculation of the processed depths in the investigated site of Altamira Cave.

Selected Processing Flow
Post-acquisition processing procedures were applied to the raw field GPR data sets to produce amplitude maps and rendered 3D images with all reflection profiles in the grids [61,62,69].In this study, a data-processing procedure was applied to the raw data sets using RADAN 7.6 software (Geophysical Survey Systems, Inc., GSSI).The main processing steps are shown in Figure 9.The first step was a 1D processing consisting of zerotime correction (time-zero adjustment).After this, the 2D processing was applied to the reflection profiles according to these steps: (i) background removal; (ii) bandpass filters; (iii) linear amplitude gains; (iv) a time-to-depth conversion based on the calculated velocity average for the overlying rock strata.Next, each processed reflector profile was aligned within the grids described, with corresponding topographic corrections, to produce 3D horizontal amplitude maps and isosurface images (isoamplitude surface rendering).In the field campaign, 454 profiles were recorded: 329 profiles totalling 20,850 m 2 were surveyed with the 100 MHz antenna, 76 profiles totalling 850 m 2 with the 400 MHz antenna, and 49 profiles totalling 550 m 2 with the 900 MHz antenna.

Data Processing Velocity Estimation
The GPR reflection data was acquired in the same week under very similar weather conditions with no rainfall.Thus, the subsurface moisture conditions in this investigated area could be considered almost constant for the investigated depth during the acquisition of the GPR field data.This velocity estimation assumption was made for the velocity-todepth conversions, considering that possible changes in subsurface conditions could cause velocity variations that would affect these velocity-to-depth conversions [61,69].
The average velocity of the GPR wave was determined by the hyperbola fitting method on a set of hyperbolas recorded in different profiles, obtaining an average velocity value of 11.17 cm/ns.The dielectric permittivity (ε) was calculated to be 7.2, according to the following equation for low-loss media [61,62]: where v is the electromagnetic wave velocity and c is the velocity of light in free space (c ≈ 29.97 cm/ns).This value was applied for the calculation of the processed depths in the investigated site of Altamira Cave.

Selected Processing Flow
Post-acquisition processing procedures were applied to the raw field GPR data sets to produce amplitude maps and rendered 3D images with all reflection profiles in the grids [61,62,69].In this study, a data-processing procedure was applied to the raw data sets using RADAN 7.6 software (Geophysical Survey Systems, Inc., GSSI).The main processing steps are shown in Figure 9.The first step was a 1D processing consisting of zero-time correction (time-zero adjustment).After this, the 2D processing was applied to the reflection profiles according to these steps: (i) background removal; (ii) bandpass filters; (iii) linear amplitude gains; (iv) a time-to-depth conversion based on the calculated velocity average for the overlying rock strata.Next, each processed reflector profile was aligned within the grids described, with corresponding topographic corrections, to produce 3D horizontal amplitude maps and isosurface images (isoamplitude surface rendering).

Integration of UAV and 3DTLS
First, the UAV terrain model was integrated with the 3DTLS model of the cave because all information is georeferenced.From this, the surface position of the galleries and the proximity of the ceiling to the exterior terrain could be analyzed (Figure 10).A cartography was created with the overlying layer thickness of the cave (Figure 11) [29].These maps consist of a raster with a resolution of 1 × 1 cm 2 , in which each pixel has the minimum distance between the ceiling of the cave and the ground recorded by the UAV flight for each pair of XY coordinates [70].

Integration of UAV and 3DTLS
First, the UAV terrain model was integrated with the 3DTLS model of the cave because all information is georeferenced.From this, the surface position of the galleries and the proximity of the ceiling to the exterior terrain could be analyzed (Figure 10).

Integration of UAV and 3DTLS
First, the UAV terrain model was integrated with the 3DTLS model of the cave because all information is georeferenced.From this, the surface position of the galleries and the proximity of the ceiling to the exterior terrain could be analyzed (Figure 10).A cartography was created with the overlying layer thickness of the cave (Figure 11) [29].These maps consist of a raster with a resolution of 1 × 1 cm 2 , in which each pixel has the minimum distance between the ceiling of the cave and the ground recorded by the UAV flight for each pair of XY coordinates [70].A cartography was created with the overlying layer thickness of the cave (Figure 11) [29].These maps consist of a raster with a resolution of 1 × 1 cm 2 , in which each pixel has the minimum distance between the ceiling of the cave and the ground recorded by the UAV flight for each pair of XY coordinates [70].

Intregration of UAV and GPR
After processing the data recorded with the three central frequency antennas (100 MHz, 400 MHz and 900 MHz), we obtained 3D isosurface models with significant discontinuities in the overlying layer of the Altamira Cave bed and with higher resolution in the layer overlying the Polychrome Hall.

Results from the Overlying Layer of the Altamira Cave Site Derived from the 100 MHz Antenna
Subvertical reflectors and vertical diffractions were mapped and recorded in the GPR profiles that correspond to internal discontinuities in the rock mass (fractures, joints, detachments, bedding planes, etc.).Some of the vertical/subvertical discontinuities, with great in-depth development, were only a few meters from the surface.These vertical/subvertical discontinuities represent important accesses/pathways for the exchange processes between the interior of the cave and the exterior.Mainly, these vertical/subvertical discontinuities were detected in the following areas (Figure 12):  The surrounding sinkholes (cave entrance, areas of Hoya Hall and El Pozo).

Intregration of UAV and GPR
After processing the data recorded with the three central frequency antennas (100 MHz, 400 MHz and 900 MHz), we obtained 3D isosurface models with significant discontinuities in the overlying layer of the Altamira Cave bed and with higher resolution in the layer overlying the Polychrome Hall.

Results from the Overlying Layer of the Altamira Cave Site Derived from the 100 MHz Antenna
Subvertical reflectors and vertical diffractions were mapped and recorded in the GPR profiles that correspond to internal discontinuities in the rock mass (fractures, joints, detachments, bedding planes, etc.).Some of the vertical/subvertical discontinuities, with great in-depth development, were only a few meters from the surface.These vertical/subvertical discontinuities represent important accesses/pathways for the exchange processes between the interior of the cave and the exterior.Mainly, these vertical/subvertical discontinuities were detected in the following areas (Figure 12):

•
The surrounding sinkholes (cave entrance, areas of Hoya Hall and El Pozo).

•
The Hoya Hall.

•
The artificial walls.

•
An area to the NW where a house and a cowshed with its dunghill were located.

•
The end zone of the westernmost branch line.
In addition, this GPR survey pinpointed several sinkholes in two areas (Figure 12):

•
To the northwest of the Altamira Cave, specifically on the road to the south of the Museum building and where buildings and constructions demolished (house, bar, stable, and a dunghill) were located.

•
In the end zone of the westernmost branch line of Altamira Cave.

•
The Hoya Hall.

•
The artificial walls.

•
An area to the NW where a house and a cowshed with its dunghill were located.

•
The end zone of the westernmost branch line.
In addition, this GPR survey pinpointed several sinkholes in two areas (Figure 12): • To the northwest of the Altamira Cave, specifically on the road to the south of the Museum building and where buildings and constructions demolished (house, bar, stable, and a dunghill) were located.

•
In the end zone of the westernmost branch line of Altamira Cave.The existence of vertical/subvertical discontinuities (fractures) detected in the Hoya Hall area would justify the construction of the artificial walls at the beginning of the 20th century.This mapping has revealed discontinuities in specific areas of the Cave de Altamira that suggest that in past times, the extension of the Cave was greater than it is today and that it had several entrances/accesses.3.2.2.Results from the Overlying Layer of the Polychrome Hall Derived from the 400 MHz and 900 MHz Antennae The 3D model from the data recorded with the 400 MHz center frequency antenna indicated two levels containing a large number of high amplitude reflectors and vertical diffractions in the first few meters of depth of the overlying layer (Figure 13):

•
A first level from the surface to a depth of about 1.20 m, which corresponds to the lapiés zone.

•
A second level located in a variable interval between 1.70 m and 2.80 m depth, which may correspond mainly to the existence of fractures, joints, and an extended zone of detachments.
Remote Sens. 2023, 15, x FOR PEER REVIEW 15 of 25 This mapping has revealed discontinuities in specific areas of the Cave de Altamira that suggest that in past times, the extension of the Cave was greater than it is today and that it had several entrances/accesses.

Results from the Overlying Layer of the Polychrome Hall Derived from the 400 MHz and 900 MHz Antennae
The 3D model from the data recorded with the 400 MHz center frequency antenna indicated two levels containing a large number of high amplitude reflectors and vertical diffractions in the first few meters of depth of the overlying layer (Figure 13):  A first level from the surface to a depth of about 1.20 m, which corresponds to the lapiés zone. A second level located in a variable interval between 1.70 m and 2.80 m depth, which may correspond mainly to the existence of fractures, joints, and an extended zone of detachments.
A presence of a number of vertical/subvertical discontinuities (fractures) connecting these two levels and extending in depth (>3.80 m depth) was also detected.These mapped discontinuities represent significant access/pathways for the exchange processes between the interior of the Polychrome Hall and its external environment.The second level was detected only to the north and east of the projection of the Polychrome Hall (Figure 14).This level could represent an area of water load located at a higher level than the ceiling of the Polychrome Hall, i.e., there is a probability of a local and temporary generation of a hanging aquifer where water from the surrounding area would build up on possible impermeable levels at its base.It is also necessary to consider the possibility that in times of high exterior-interior exchange, many particles deposited there may end up inside the Polychrome Hall.A presence of a number of vertical/subvertical discontinuities (fractures) connecting these two levels and extending in depth (>3.80 m depth) was also detected.These mapped discontinuities represent significant access/pathways for the exchange processes between the interior of the Polychrome Hall and its external environment.
The second level was detected only to the north and east of the projection of the Polychrome Hall (Figure 14).This level could represent an area of water load located at a higher level than the ceiling of the Polychrome Hall, i.e., there is a probability of a local and temporary generation of a hanging aquifer where water from the surrounding area would build up on possible impermeable levels at its base.It is also necessary to consider the possibility that in times of high exterior-interior exchange, many particles deposited there may end up inside the Polychrome Hall.In addition, this second level could contribute to the direct and lateral infiltration (impluvium) of water into the horizons that make up the geological structure of the overlying layer of Polychrome Hall.This is the case in three of the vertical fracture zones observed and coincides with the temporary contributions of water in the walls of the Polychrome Hall, with two of them being located in its north wall and one of them in its south wall, as shown in Figure 15.In addition, this second level could contribute to the direct and lateral infiltration (impluvium) of water into the horizons that make up the geological structure of the overlying layer of Polychrome Hall.This is the case in three of the vertical fracture zones observed and coincides with the temporary contributions of water in the walls of the Polychrome Hall, with two of them being located in its north wall and one of them in its south wall, as shown in Figure 15.In addition, the existence of this detachment level (level of high-amplitude reflectors), which borders the Polychrome Hall to the north and east, could be related to

•
the decision at the beginning of the 20th century to construct the artificial walls to the north and east of the Polychrome Hall, or • the current topography of the surface above Polychrome Hall, reflecting a change in orientation and slope that tilts towards the NE on the surface, which coincides with the surface projection of this detected level of high amplitude reflectors (Figure 16), or  In summary, we observed that the overlayer of Polychrome Hall is more affected by fracturing in its northern and eastern areas.Similarly, its fracturing increases from west to east.
The 900 MHz center frequency antenna was used for the 3D study of the lapiés area of the Polychrome Hall.In this study, a large number of high-amplitude reflectors were recorded and mapped from the surface to about 1.20 m depth.These high-amplitude reflectors correspond to the rugged water-modeled limestone surface, consisting of rock pinnacles separated by grooves.The lineaments of this groove complex in the lapiés derived from 3D mapping agree with the predominant directions of the joint system (N45, N68, and N168) mapped in the Polychrome Hall rock outcrops [68].(Figure 17).This lapiés area represents a permeable soil system that allows both runoff and infiltration of water, as well as the flow of gases.In summary, we observed that the overlayer of Polychrome Hall is more affected by fracturing in its northern and eastern areas.Similarly, its fracturing increases from west to east.
The 900 MHz center frequency antenna was used for the 3D study of the lapiés area of the Polychrome Hall.In this study, a large number of high-amplitude reflectors were recorded and mapped from the surface to about 1.20 m depth.These high-amplitude reflectors correspond to the rugged water-modeled limestone surface, consisting of rock pinnacles separated by grooves.The lineaments of this groove complex in the lapiés derived from 3D mapping agree with the predominant directions of the joint system (N45, N68, and N168) mapped in the Polychrome Hall rock outcrops [68].(Figure 17).This lapiés area represents a permeable soil system that allows both runoff and infiltration of water, as well as the flow of gases.

Visualization of the Integration of UAV, 3DTLS, and GPR
The 3D GPR model derived from the 100 MHz antenna recordings provides information about the most important discontinuities in the Altamira Cave karst system up to a depth of about 8 m (Figure 18).The 3D GPR model derived from the 100 MHz antenna recordings provides information about the most important discontinuities in the Altamira Cave karst system up to a depth of about 8 m (Figure 18).This model, showing the 3D amplitude isosurfaces, has been integrated with the UAV and 3DTLS models in Blender, an open-source 3D computer graphics software [72], to manage data and determine the spatial relationships between them by measuring or making sections (Figure 19) with variable thickness.This model, showing the 3D amplitude isosurfaces, has been integrated with the UAV and 3DTLS models in Blender, an open-source 3D computer graphics software [72], to manage data and determine the spatial relationships between them by measuring or making sections (Figure 19) with variable thickness.
A video (Figure 20) showing the remarkable discontinuities recorded during the geophysical campaign in the karstic formation of the Altamira Cave was also rendered.The 3D GPR model derived from the 100 MHz antenna recordings provides information about the most important discontinuities in the Altamira Cave karst system up to a depth of about 8 m (Figure 18).This model, showing the 3D amplitude isosurfaces, has been integrated with the UAV and 3DTLS models in Blender, an open-source 3D computer graphics software [72], to manage data and determine the spatial relationships between them by measuring or making sections (Figure 19) with variable thickness.A video (Figure 20) showing the remarkable discontinuities recorded during the geophysical campaign in the karstic formation of the Altamira Cave was also rendered.

Discussion
Geologically speaking, Altamira Cave is in a senile state with a clear tendency to disappear through natural evolution, given that the processes of destruction are more important than those of sedimentation [26].The existence of an important network of fissures, fractures, joints, and displacements present from the top of the cave to the cave ceilings themselves affects cave stability and conditions important changes in the circulation of infiltration water, generating relevant alterations on the painting such as flaking, washing/erosion, and the concretion of carbonates [73,74].
The criteria used for the management and custody of Altamira Cave have led to the application of measures aimed at protecting the cave, its paintings, and its surroundings to guarantee their correct conservation.The research and conservation work in the cave, within the framework of its Preventive Conservation Plan, is aimed at gaining a better understanding of the active deterioration processes to establish a detailed diagnostic procedure that facilitates indirect conservation actions.
Accurate georeferencing is key for integrating remote-sensing techniques in order to derive useful information from the cultural heritage, since it allows us to see the relationships between different parts of the cave with the exterior and how the soil is preserved.

Discussion
Geologically speaking, Altamira Cave is in a senile state with a clear tendency to disappear through natural evolution, given that the processes of destruction are more important than those of sedimentation [26].The existence of an important network of fissures, fractures, joints, and displacements present from the top of the cave to the cave ceilings themselves affects cave stability and conditions important changes in the circulation of infiltration water, generating relevant alterations on the painting such as flaking, washing/erosion, and the concretion of carbonates [73,74].
The criteria used for the management and custody of Altamira Cave have led to the application of measures aimed at protecting the cave, its paintings, and its surroundings to guarantee their correct conservation.The research and conservation work in the cave, within the framework of its Preventive Conservation Plan, is aimed at gaining a better understanding of the active deterioration processes to establish a detailed diagnostic procedure that facilitates indirect conservation actions.
Accurate georeferencing is key for integrating remote-sensing techniques in order to derive useful information from the cultural heritage, since it allows us to see the relationships between different parts of the cave with the exterior and how the soil is preserved.The study carried out on the mountain where Altamira Cave itself is located posed a research challenge of great relevance for the cave's conservation: to determine the exact access routes and the trajectories of the infiltration water which, through the lines of fractures, fissures, and joints, access the different parts of the cave, especially the Polychrome Hall.
Integrating GNSS, 3DTLS, UAV photogrammetry and GPR permits a rapid, accurate, and reliable recording of complex features such as caves or cavities [18].This information can be used to derive cartography such as floor plans, contour lines, longitudinal and cross sections, and three-dimensional analysis such as the height or thickness of galleries, or analyze and map the system of fractures from the exterior to almost the cave ceiling.Very significant vertical and subvertical fractures were detected, as well as the connections between these fractures and joints, establishing the access routes for infiltration water and relating them to the drip points located on the ceilings of the cave.All this information on moisture that is directly or indirectly associated with the main fractures and fissures mapped has been linked to many of the drip points involved in processes of erosion, flaking, disintegration, and washing away by dissolution of some of the mineral components present in the cave ceilings [75].
This model integration supports the decision-making of cave managers by enabling them to engage with advanced predictive models and simulations within a knowledgebased system.

Figure 1 .
Figure 1.Location map and the surroundings of Altamira Cave.The cave is located in Santillana del Mar, Cantabria, Spain.

Figure 1 .
Figure 1.Location map and the surroundings of Altamira Cave.The cave is located in Santillana del Mar, Cantabria, Spain.

Figure 2 .
Figure 2. (a) Map of the interior of the Cave including gallery names, artificial walls, pillars, and stairs.(b) Map of the traverse and reference checkerboards, including the traverse station from which they radiate.

Figure 3 .
Figure 3. General workflow scheme followed in this study.

Figure 3 .
Figure 3. General workflow scheme followed in this study.

Figure 4 .
Figure 4. View of the 3D model derived from UAV photogrammetric flight of the surroundings of Altamira Cave.

Figure 4 .
Figure 4. View of the 3D model derived from UAV photogrammetric flight of the surroundings of Altamira Cave.

Figure 5 .
Figure 5.The location of GPR profile grid (orthogonal 2D profiles spaced every 3 m) in the Altamira Cave site, carried out with the 100 MHz antenna.

Figure 5 .
Figure 5.The location of GPR profile grid (orthogonal 2D profiles spaced every 3 m) in the Altamira Cave site, carried out with the 100 MHz antenna.

Figure 6 .
Figure 6.Radargram recorded with the 100 MHz antenna.Profile Z2-85 above shows discontinuities in the rock mass: several sinkholes (indicated by white ovals) and vertical fractures.This radargram also shows the attenuation of the electromagnetic wave in the last meters of depth.

Figure 6 .
Figure 6.Radargram recorded with the 100 MHz antenna.Profile Z2-85 above shows discontinuities in the rock mass: several sinkholes (indicated by white ovals) and vertical fractures.This radargram also shows the attenuation of the electromagnetic wave in the last meters of depth.

Figure 7 .Figure 7 .Figure 7 .Figure 8 .
Figure 7.The location of GPR profile grids (orthogonal 2D profiles spaced every 1 m) in the overlying rock strata and the lapiés zone of the Polychrome Hall ceiling, carried out with the (a) 400 MHz and (b) 900 MHz antennas.

Figure 8 .
Figure 8.(a) Location and (b) radargram (profile T-97) recorded with the 400 MHz antenna, in which two levels of discontinuities in the rock mass were detected.The first and more superficial level corresponds to the lapiés (delimited with a black dashed line), and the second corresponds to a zone of fractures, joints, and detachments (delimited with a yellow dashed line).Vertical discontinuities corresponding to vertical fractures were also detected.(c) Radargram (profile T-97) recorded with the 900 MHz antenna, in which discontinuities with better resolution were detected at the lapiés level (delimited with a black dashed line).Vertical discontinuities corresponding to vertical fractures and rock pinnacles separated by grooves were also detected.

Figure 8 .
Figure 8.(a) Location and (b) radargram (profile T-97) recorded with the 400 MHz antenna, in which two levels of discontinuities in the rock mass were detected.The first and more superficial level corresponds to the lapiés (delimited with a black dashed line), and the second corresponds to a zone of fractures, joints, and detachments (delimited with a yellow dashed line).Vertical discontinuities corresponding to vertical fractures were also detected.(c) Radargram (profile T-97) recorded with the 900 MHz antenna, in which discontinuities with better resolution were detected at the lapiés level (delimited with a black dashed line).Vertical discontinuities corresponding to vertical fractures and rock pinnacles separated by grooves were also detected.

Figure 9 .
Figure 9. Flowchart of the main processing steps applied to raw reflection data collected by GPR in this study.

Figure 10 .
Figure 10.Visualization of the 3D model of the cave below a low-detail mesh obtained from UAV model.

Figure 9 .
Figure 9. Flowchart of the main processing steps applied to raw reflection data collected by GPR in this study.

25 Figure 9 .
Figure 9. Flowchart of the main processing steps applied to raw reflection data collected by GPR in this study.

Figure 10 .
Figure 10.Visualization of the 3D model of the cave below a low-detail mesh obtained from UAV model.

Figure 10 .
Figure 10.Visualization of the 3D model of the cave below a low-detail mesh obtained from UAV model.

Figure 12 .
Figure 12.Zonal distribution map of subvertical reflectors and vertical diffractions (fractures) from GPR profiles in the rock massif of Altamira Cave: (a) the surrounding sinkholes, cave entrance, areas of Hoya Hall and Pozo Hall; (b) the Polychrome Hall; (c) the Cave of the Stalactites; (d) the Hoya Hall; (e) the artificial walls; (f) an area to the NW where a house and a cowshed with its dunghill were located; (g) the end zone of the westernmost branch line.Sinkhole areas pinpointed by GPR in the rock massif of Altamira Cave are indicated by yellow ovals.

Figure 12 .
Figure 12.Zonal distribution map of subvertical reflectors and vertical diffractions (fractures) from GPR profiles in the rock massif of Altamira Cave: (a) the surrounding sinkholes, cave entrance, areas of Hoya Hall and Pozo Hall; (b) the Polychrome Hall; (c) the Cave of the Stalactites; (d) the Hoya Hall; (e) the artificial walls; (f) an area to the NW where a house and a cowshed with its dunghill were located; (g) the end zone of the westernmost branch line.Sinkhole areas pinpointed by GPR in the rock massif of Altamira Cave are indicated by yellow ovals.The existence of vertical/subvertical discontinuities (fractures) detected in the Hoya Hall area would justify the construction of the artificial walls at the beginning of the 20th century.

Figure 13 .
Figure 13.3D amplitude isosurface model of the main discontinuities in the overlying layer of the Polychrome Hall with the 400 MHz antenna, showing the two levels of discontinuities (lapiés zone-fractures, joints, and detachment zone) and the vertical fractures.

Figure 13 .
Figure 13.3D amplitude isosurface model of the main discontinuities in the overlying layer of the Polychrome Hall with the 400 MHz antenna, showing the two levels of discontinuities (lapiés zone-fractures, joints, and detachment zone) and the vertical fractures.

Figure 14 .
Figure 14.Amplitude slice, 2.20 m depth with the 400 MHz antenna, showing the fracture, joint, and detachment zone (marked in white) on the vertical projection of the Polychrome Hall.This zone of discontinuities is located only to the north and east of the Polychrome Hall projection.

Figure 14 .
Figure 14.Amplitude slice, 2.20 m depth with the 400 MHz antenna, showing the fracture, joint, and detachment zone (marked in white) on the vertical projection of the Polychrome Hall.This zone of discontinuities is located only to the north and east of the Polychrome Hall projection.

Figure 15 .
Figure 15.Amplitude slice, 3.50 m depth with the 400 MHz antenna, showing with white arrows the locations of the projections of the vertical/subvertical fractures that correspond to the 3 zones of temporary water inflow observed in the interior of the Polychrome Hall.Two of them are on the north wall and one is on the south wall.In addition, the existence of this detachment level (level of high-amplitude reflectors), which borders the Polychrome Hall to the north and east, could be related to o the decision at the beginning of the 20th century to construct the artificial walls to the north and east of the Polychrome Hall, or o the current topography of the surface above Polychrome Hall, reflecting a change in orientation and slope that tilts towards the NE on the surface, which coincides with the surface projection of this detected level of high amplitude reflectors (Figure 16), or

Figure 15 .
Figure 15.Amplitude slice, 3.50 m depth with the 400 MHz antenna, showing with white arrows the locations of the projections of the vertical/subvertical fractures that correspond to the 3 zones of temporary water inflow observed in the interior of the Polychrome Hall.Two of them are on the north wall and one is on the south wall.

•Figure 16 .
Figure 16.Projection of the Polychrome Hall on the topographic map, showing a change of orientation and slope (indicated with black arrow) that tilts to the NE (a), and which also corresponds to the projection of the second level of discontinuities in its amplitude depth slice (2.20 m depth) (b).

Figure 16 .
Figure 16.Projection of the Polychrome Hall on the topographic map, showing a change of orientation and slope (indicated with black arrow) that tilts to the NE (a), and which also corresponds to the projection of the second level of discontinuities in its amplitude depth slice (2.20 m depth) (b).

Figure 17 .
Figure 17.Amplitude slice, 0.45 m depth with the 900 MHz antenna, showing the mapped lapiés groove system in the overlying layer of the Polychrome Hall.It can be seen that the lineaments of this system of lapiés grooves, derived from the 3D mapping, agree with the orientations of the main lineaments of the joints (N45, N68, and N168)[67].

Figure 17 .
Figure 17.Amplitude slice, 0.45 m depth with the 900 MHz antenna, showing the mapped lapiés groove system in the overlying layer of the Polychrome Hall.It can be seen that the lineaments of this system of lapiés grooves, derived from the 3D mapping, agree with the orientations of the main lineaments of the joints (N45, N68, and N168)[67].

Figure 18 .
Figure 18.Image of integrating the UAV model, 3DTLS and GPR (100 MHz antenna), showing the discontinuities in the Polychrome overlying layer (in light brown) and the backside of the Polychrome Hall (in dark brown).

Figure 18 .
Figure 18.Image of integrating the UAV model, 3DTLS and GPR (100 MHz antenna), showing the discontinuities in the Polychrome overlying layer (in light brown) and the backside of the Polychrome Hall (in dark brown).

Figure 18 .
Figure 18.Image of integrating the UAV model, 3DTLS and GPR (100 MHz antenna), showing the discontinuities in the Polychrome overlying layer (in light brown) and the backside of the Polychrome Hall (in dark brown).

Figure 19 . 25 Figure 19 .
Figure 19.(a) Section of the integrated model in the Entrance Hall (1) and Polychrome Hall (2) area showing terrain, GPR (brown), and 3DTLS (grey) models.(b) Section with 6 m far plane of the integrated model in the Entrance Hall (1) and Polychrome Hall (2) area showing terrain, GPR (brown), and 3DTLS (grey) models.