Extreme Sea Surges, Tsunamis and Pluvial Flooding Events during the Last ~1000 Years in the Semi-Arid Wetland, Coquimbo Chile

Abstract

The coast of Chile has been exposed to marine submersion events from storm surges, tsunamis and flooding due to heavy rains. We present evidence of these events using sedimentary records that cover the last 1000 years in the Pachingo wetland. Two sediment cores were analyzed for granulometry, XRF, pollen, diatoms and TOC. Three extreme events produced by marine submersion and three by pluvial flooding during El Niño episodes were identified. Geochronology was determined using a conventional dating method using ¹⁴C, ²¹⁰Pbxs and ¹³⁷Cs). The older marine event (E1) was heavier, identified by a coarser grain size, high content of seashells, greater amount of gravel and the presence of two rip-up clasts, which seems to fit with the tsunami of 1420 Cal AD. The other two events (E3 and E5) may correspond to the 1922 (E3) tsunami and the 1984 (E5) storm waves, corroborated with a nearshore wave simulation model for this period (SWAM). On the other hand, the three flood events (E2, E4, E6) all occurred during episodes of El Niño in 1997 (E6), 1957 (E4) and 1600 (E6), represented by layers of fine-grain sands and wood charcoal remains.

1. Introduction

Chile has been frequently affected by extreme telluric events such as large-scale earthquakes and tsunamis or extreme meteorological events such as storms and heavy rains. The latest significant extreme telluric events are (i) three more significant earthquakes: in Valdivia in 1960, with a magnitude of 9.5° Mw; (ii) in El Maule in 2010, with a magnitude of 8.8° Mw; and (iii) in Illapel in 2015, with a magnitude of 8.3° Mw. In all cases, a tsunami occurred a few minutes later [1,2,3,4,5,6]. The one that has caused the most human and economic losses is the 1960 tsunami, in Chile alone, it caused the death of 1700 people, 3000 wounded, 2,000,000 homeless victims. reflecting in an economic loss of 550 million dollars (CIGIDEN: Centro de Investigación para la Gestión Integrada del Riesgo de Desastres/Chile). In the 2015 tsunami, in Coquimbo, deposits with a thickness that varies between 10 and 50 cm were observed [7]. The latest extreme meteorological events were (i) 21–22 May 1957, with rainfall of 77.4 mm per month affecting La Serena city. There were overflows of the Elqui river, and more than 200 people were affected by the flood; the main routes were cuts, and two people died. (ii) 11 June 1997, rainfall of 104.9 mm per month. This time the rain was accompanied by hail, thunderstorms and winds of 90–100 km per hour. It caused the activation of two streams that devastated eight houses and caused the death of two people (/https://explorador.cr2.cl/ (access date: 20 October 2021)). (iii) 22 January 2015, which produced mudslides affecting the three most extensive areas of northern Chile in the pre-mountain range sectors, leaving half of Copiapo city (Atacama region) buried under ~2 m of mud, with a high human and economic cost of $1500 million dollars (MINVU-Chile [7,8,9].

The latest tsunami events to have affected the Coquimbo region occurred in 1730, 1880, 1922, 1943 and 2015, after large ruptures of the fault line [2,10,11,12,13,14,15,16,17,18,19,20]. Of these, the most destructive according to historical records occurred in 1922, which affected 400 km of coastline from the Atacama to the Coquimbo regions (27–30° S) [16,20,21,22] and caused extensive damage in an area stretching approximately from Caldera to Coquimbo (27–30° S). More than 1000 deaths resulted from the earthquake, and the tsunami killed hundreds of people in coastal areas, mostly in Coquimbo. The range of total damage was estimated to be 5–25 million USD [23]. Other marine events affecting the coastal area have been the great wave storms. The main events occurred on 10 August 1965, 10 July 1984, and 8 August 2015, causing damaged coastal infrastructure and eroding beaches and dunes [24,25,26]. They were categorized as M4 to M5, according to the official wave intensity scale (https://marejadas.uv.cl/index.php/categorias/folleto-categorias (access date: 23 July 2021)), or 6 to 7 on the scale proposed by Campos-Caba (2016) (https://oleaje.uv.cl/marejadas.html (access date: 5 October 2021)). These categories indicate “the evacuation of the coast is imminent, the structures are severely damaged or destroyed, the persistent overpass generates flows in walks and streets, causing significant damage or destruction of properties”.

Therefore, Coquimbo has experienced numerous extreme coastal events, caused by meteorological conditions as well as tsunami events. These have caused casualties and economic damages [27,28], and thus it is necessary to study the frequency and intensity of past events to predict future trends, determine their frequency and prevent risks. However, historical records and meteorological data limit the analysis to only a few centuries [29,30]. Rainfall is recorded by the General Directorate of Water (DGA; https://explorador.cr2.cl/, access date: 20 October 2021), and for the Pachingo wetland, the data have been collected since 1930. In the case of storms, wave record data have been collected since 1823 (Campos-Caba et al., 2016/https://oleaje.uv.cl/marejadas.html, access date: 20 October 2021; https://marejadas.uv.cl/, access date: 20 October 2021) and tsunami records for some centuries, since 1562, although the oldest records are based only on historical records in habited areas (https://www.proteccioncivil.es/catalogo/naturales/jornada-maremotos/documentacion/docu2.pdf/, access date: 20 October 2021; CERESIS: Regional Seismology Center for South America. Catalogue of Earthquakes for South America. 1985; NOAA/https://www.ngdc.noaa.gov/hazel/view/hazards/tsunami/event-data?maxYear=2015&minYear=1200&country=CHILE, access date: 20 October 2021). Tsunami records begin after the Spanish invasion in 1535 [20]. Meanwhile, geological records contribute to the reconstitution of climatic weather events, intense storms and tsunami activity at a longer scale [30]. In this sense, there are other environments with similar suitability for such studies in coastal areas are wetlands with sandy barriers. These natural environments have been highly affected by changes in water inputs controlled by the amount of precipitation and by the geo-morphological alterations of the basin and coast, caused by the sea-level changes, tidal waves, tectonic movements, tsunamis and flooding events [31,32,33,34,35]. These changes can be observed in the sediment stratification because there is scarce post-depositional remobilization due to the roots’ mechanical action, which retains the organic and clayey particulate material [36,37]. The interpretation of proxies in these sedimentary environments is relevant to establishing meteorological events and extreme sea events of the past [31,38,39]. The identification of different depositional environments formed by the recurrence of overwash layers and associated contents allows the reconstruction of marine paleo-events [31,32,35,39,40,41,42,43,44,45,46,47,48,49,50,51,52]. In this study, we reconstruct past extreme sea events (storm waves or tsunamis) and pluvial flooding events (extreme rains) using a multiproxy analysis based on sedimentological, biological and geochemical data in the Pachingo wetlands. In these sediments, we analyzed wood charcoal and carbonates for the age estimations. In addition, a swam model was used to understand how storms waves affect the coastal area of the Pachingo wetlands. Results were then compared with the historical extreme sea event records available at Pachingo to identify the events that may have reached or submerged the sites.

2. Study Area

Pachingo’s wetland (30°18′36.00″ S; 71°34′17.21″ W) is near to Tongoy Bay in the Coquimbo region, Chile. It is a glacio-eustatic bay formed by the abrasive action of the successive interglacial transgressions that occurred during the Quaternary [53,54]. These movements formed the stepped littoral marine terraces, caused by coastal elevation, where three wetlands were formed—Salinas Chicas, Salinas Grande and Pachingo—with a progradation rate of ca. 0.14 m [14,32,55]. The Coquimbo region is a semi-arid zone, strongly affected by the warm phase of El Niño Southern Oscillation (ENSO) when rainfall increases strongly compared to neutral periods or during the cold phase of La Niña. It causes streams to grow that generally are dry or have low flows [56,57,58]. The intense flood events recorded by the Meteorological Directorate of Chile are from 1957, 1965, 1972, 1987, 1991, 1992, 2015 and 2017, all of which were related to the El niño phenomenon [56]. The wetland has a permanent shallow stream and is separated from the sea by a dune barrier [55,59,60], which disappears in cases of heavy rains that generate large floods, being temporarily connected to the sea [61,62]. The creeks are fed by the runoff from small river inputs (Elqui and Limari), groundwater flow and precipitation in the basin [61,62,63]. Studies in different terrace levels of the Pachingo’s ravine [55] showed evidence of the occurrence of torrential rivers in several periods during the mid-to-late Holocene, and repeated tectonic elevation of the coastal segment occurred during the Holocene [32].

The wetland’s total area is ~0.3 km², and it is at a height of 0.5 to 1.0 MSL (mean sea level). The climate is characterized by high humidity (85%), the temperature varies between 11–22 °C, and the annual temperature average is 14.7 °C. Precipitation is concentrated in the winter months (May to August), and it does not exceed 75 mm per year (annual average last 10 years: 90 a 100 mm) except during El Niño (200 to 400 mm) [56,57,58,59,64].

3. Methodology

Sediment cores PT1 (50 cm long) and PT4 (47 cm long), were obtained in 2014 in the Pachingo wetland at ~850 m and 920 m from the coastline, respectively—an area currently flooded by tsunamis [1] (Figure 1). The core was obtained with a manual core under 20 and 50 cm of water depth, respectively. X-rays were taken for the PT1 core in the San Juan de Dios Hospital of La Serena, Chile, performed with the methodology used for femur and spinal column inspection. Then, the cores were sectioned every 1 cm, and the samples were stored in sterilized Whirlpak^® bags, frozen and subsequently freeze-dried.

Visual stratigraphic examination was conducted for both cores. The PT1 core was analyzed for sedimentology (grain size, Kurtosis and Skewness, magnetic susceptibility), geochemistry (elemental composition, TOC, biogenic opal, ¹⁵C,¹³N) and biological proxies (Diatoms and pollen). Diatoms and pollen were identified, and their abundance was estimated. The Geochronology analyses ¹³⁷Cs, ²¹⁰Pb and ¹⁴C were used. The core PT4 was only analyzed for granulometry and geochemistry to compare the events identified at both cores.

3.1. Lithology and Stratigraphy

The lithological units were visually characterized according to the Munsell chart scale, X-rays, granulometry and the types of contact between units. These units were represented with a schematic log. Magnetic susceptibility (SI × 10⁻⁸) was measured with a Bartington Susceptibility Meter MS2E in the Sedimentology Laboratory at Centro Eula, Universidad de Concepción. Samples were measured three times, and the values were expressed as mean values.

The grain size was determined using a Beckman–Coulter LS13320 laser diffraction particle size analyzer (Géosciences Montpellier Laboratory, Montpellier, France). The analysis was carried out in each centimeter, considering particles <2 mm in diameter using a sieve. The particle sizes were classified according to the scale of Folk and Ward (1957).

Sorting, skewness and kurtosis were evaluated using the GRADISTAT statistical software using the method moments logarithmic according to Blott and Pye, (2001); which includes all particle size spectra.

In addition to laser particle size analysis (<2 mm), we performed sieve analyses (0.63, 0.08, 0.1, 0.125, 0.16, 0.2, 0.25, 0.315, 0.4, 0.5, 0.63, 0.8, 1, 1.25, 1.6, 2, 2.5, 3.15, 4, 5 and 6.3 mm) to specify more specifically the distribution of particles greater than 2 mm on two key samples (E1 and E5). However, no statistical analysis was performed on this coarser fraction.

3.2. Geochemical Characterization

Elemental analyses were performed by X-ray fluorescence (XRF) spectrometry for each centimeter on the surface of the PT1 core using a hand-held Niton XL3t spectrometer (pXRF; CRAHAM laboratory, Université de Caen-Normandie, Caen, France). The analysis considered Ba, K, Rb, Sr, Ca, P, Mg, Mn, Zr, Fe, Zn, Ti, Ni and V. Principal Components Analysis (PCA) was used to identify correlated variables between geochemical poles using the XLStat Program and to identify the most representative element of each geochemical family.

3.3. Geochronology

Radiochronological data for the last century were obtained using ²¹⁰Pbxs (activity “in excess” or unsupported) on 22 selected samples along the PT1 core. This method consists of evaluating ²¹⁰Po, in secular equilibrium with ²¹⁰Pb. The ²¹⁰Po activities were quantified by alpha spectroscopy using a CANBERRA QUAD alpha spectrometer (model 7404) until obtaining an adequate statistical count (~4–10%, error). The chemical procedure considered the digestion of ~0.5 g of sediment with concentrated acids in the presence of a chemical tracer (²⁰⁹Po, 2.22 dpm g⁻¹) and the plating on silver discs 99.9% pure for 3 h at ~75 °C according to Flynn (1968). The alpha counting was performed at the Oceanography Department of Universidad de Concepción, Chile.

Geochronology was then established based on excess activities following the CRS method described by McCaffrey and Thomson (1980); details are described in Muñoz and Salamanca (2003). Additionally, ¹³⁷Cs measurements were quantified with a BEGe 3825 gamma spectrometer at Géosciences Montpellier Laboratory (Montpellier, France).

Radiocarbon measurements were performed to obtain older ages (>200 years) on five estuarine snail shells (Heleobia sp.) and one sample of wood charcoal obtained from PT1 core. The samples were dated using conventional 14C-AMS measurements submitted to the DirectAMS (Radiocarbon Dating Services; directams.com, access date: 20 October 2021). The dates obtained were calibrated using the marine 20.14C calibration curve for shell samples and Shcal20.14C curve for southern hemisphere 14C dates for wood charcoal [65,66]. The dates obtained from shells were adjusted with the three reservoir ages calculated for the area (

Table 1. Radiocarbon details and resulted ages on shells (Heleobia sp.) and wood charcoal from the PT1 core. The dates obtained for the snail shells were calibrated using Calib 8.20, and the marine 20.14C curve, testing three reservoir ages (a, b and c); for wood charcoal age, it was calibrated using the Shcal20.14C curve [65,66].

				Fraction Modern		14C Age		Age Calib
ID	Material	Section	Lab Code	pMC	1σ Error	BP	1σ Error	BC/AD	BC/AD	BC/AD	BC/AD
Core		cm	DirectAMS					210Pb	ΔR 625 ± 46 a	ΔR 165 ± 107 b	ΔR 442 c
PT1	Shell	1	D-AMS 009072	94.27	0.29	474	25	2010	invalid age	invalid age	invalid age
PT1	Shell	8	D-AMS 009060	87.62	0.26	1062 *	24	1986	invalid age	1468–1700	1832–1950
PT1	Shell	23	D-AMS 009061	76.14	0.22	2190 *	23	-	901–1072	403–658	729–871
PT1	Wood charcoal	24	D-AMS 013087	96.68	0.34	271	28	-	1640–1671	1640–1671	1640–1671
PT1	Shell	36	D-AMS 009063	74.59	0.27	2355 *	30	-	726–902	199–491	580–700
PT1	Shell	49	D-AMS 009064	80.56	0.28	1736 *	28	-	1336–1468	858–1129	1205–1324

* reverse ages: observed in extreme event deposits; ^a Ortega et al., 2019; ^b Carré et al., 2016; ^c Muñoz et al., 2020.

).

3.4. Biogenic Components in the Sedimentary Records

Diatoms and pollen are relevant tools to reconstruct the environmental and meteorological variability of the past [67,68,69,70,71,72,73,74,75]. For the determination of diatoms, 500 mg of freeze-dried sediments were used. The samples were oxidized and mounted following the methodology of Rebolledo et al. (2005). The clean samples were observed under a Carl Zeiss microscope at 400 and 1000× using a contrast phase. Fragments containing more than half a valve were included in the count. Diatoms were identified to the lowest taxonomic level and were grouped by fresh water, brackish and marine species following [76,77,78]. The counts were expressed as diatom valves g⁻¹. Biogenic opal was determined as SiOPAL (%) according to Mortlock and Froelich (1989), and the values were estimated with a precision of ±0.5%. All laboratory and microscopy work was performed in the Paleoceanography Laboratory at the University of Concepción.

The samples for palynological determination were carried out according to the standard protocols for terrestrial sedimentary samples, including a treatment with 10% KOH solution and the dissolution of the carbonates with 5% HC. The silica-clastic fraction was dissolved with a concentrated HF solution, and the cellulose was removed with an Acetolysis reaction [79]. The samples were mounted with glycerol and permanently sealed with paraffin wax. These were analyzed using a light microscope at 40× magnification. The identification of the palynological samples was performed with the assistance of the reference samples housed in the Paleoecology laboratory of the Center for Advanced Studies in Arid Zones (CEAZA) and based on the Catalog published by Heusser (1971). The pollen count considered between 10 and 300 terrestrial grains, depending on the abundance of pollen in each sample. The counts were expressed as grains g-1. Abundant changes of each taxon were expressed in percentage values. The pollen percentage from aquatic species and ferns was calculated for the total of the terrestrial taxa. The percentage diagrams were developed with the Tilia version 2.04 software (E. Grimm, Illinois State Museum, Springfield, IL, USA).

3.5. SWAN Model

The open-source numerical model SWAN (Simulating Waves Nearshore), developed by the Delft University of Technology, was used to propagate the waves from deep waters and understand the impact of waves on the coastline. SWAN is a third-generation phase-averaged numerical wave model that simulates the effects of refraction, diffraction, break due to slenderness and bottom influence, wind effects, wave–wave interaction and white capping in coastal regions and inland waters [80]. The model solves the spectral action balance equation, allowing the calculation of the wave spectrum evolution in a time/space domain from which statistical parameters of a sea state can be determined [81].

The wave conditions of propagation were evaluated considering the peak values of the significant wave height (Hs) in deep waters, with the associated peak period (Tp) and peak direction (Dp). These variables were considered for the evaluation of events identified on 5 and 10 July 1984 (Figure 2). The parameters in deep waters are specified in

Table 2. Wave parameters in deep water for SWAN model propagation. Hs: wave height, Tp: peak period, Dp: peak direction.

Date	Hs [m]	Tp [s]	Dp [°]
5 July 1984	5.67	14.8	232.5
10 July 1984	4.87	9.2	322.5

. The statistics of waves in deep waters were obtained from the database of the project “Un Atlas de Oleaje para Chile” [82]. The deep-water wave statistics used were generated with the Wavewatch III v.4.18 model in the high-performance cluster of the Center for Mathematical Modeling of the University of Valparaíso (CIMFAV).

Regarding the configuration of the SWAN model, three structured calculation meshes were used with resolution elements of 300 × 300 (m), 100 × 100 (m) and 20 × 20 (m), reducing in size toward the coast to obtain greater precision in the area of interest. Bottom friction and dissipation due to breakage were not considered in the modeling to maintain a conservative criterion for estimating wave heights at the site of interest.

4. Results

4.1. Lithostratigraphic Description

The core PT1 has six units presented with a schematic log named from E1 to E6 (Figure 3). The deposits exhibited significant differences in texture, color and composition; the presence of shell remains and very coarse sand in the base of the core (identified as E1) was evident compared to the other units. In general, the thicker sand deposits were identified in units E1, E3 and E5 from medium sand to very coarse sand. According to the Munsell color scale, these layers varied from very dark gray (2.5Y 3/1) to dark gray (2.5Y 4/1). These three deposits (E1, E3 and E5), are generally nominated by a high density (>1.5 g/cm³), a decreasing trend of skewness (Coarse skewed <0.3), kurtosis (Platykurtic <2) and magnetic susceptibility. The three deposits present an erosive contact surface with the deposits that delimit them (Figure 4A,

Table 3. Sedimentological characteristics of meteorological event deposits along the PT1 and PT4 cores.

ID	ID	Section	Sorting	Contact		Munsell Color Scale	Other Characteristics
Core	Unit	cm		Upper	Basal
PT1	E1	30–50	Moderately sorted Very poorly sorted *	Sharp	-	2.5Y 4/1	High content of gravel, angular fragments of seashell, seashell and two rip-up clasts (30 and 40 cm)
PT1	E2	20–22	Poorly sorted	Relatively sharp	Sharp	2.5Y 2.5/1	High content of roots, wood charcoal and estuarine shell
PT1	E3	16–19	Well sorted	Relatively sharp	Relatively sharp	2.5Y 4/1	High content of gravel and fragments of angular seashell
PT1	E4	10–12	Poorly sorted	Sharp	Slightly sharp	2.5Y 2.5/1	High content of roots, round gravel, wood charcoal and estuarine shell
PT1	E5	4–8	Moderately well sorted	Relatively sharp	Sharp	2.5Y 4/1	High content of gravel and seashell fragments
PT1	E6	2–3	Moderately sorted	gradational	Relatively sharp	2.5Y 3/1	High content of roots and estuarine shell
PT4	E1	22–47	Poorly sorted	Sharp	-	2.5Y 4/1	High content of gravel, angular fragments of seashell and seashell
PT4	E2	17–20	Very poorly sorted	Relatively sharp	Sharp	2.5Y 2.5/1	High content of roots, wood charcoal and estuarine shell
PT4	E4	6–8	Poorly sorted	Slightly sharp	Slightly sharp	2.5Y 2.5/1	High content of roots and estuarine shell
PT4	E5	3–5	Poorly sorted	Relatively sharp	Sharp	2.5Y 4/1	High content of gravel, angular fragments of seashell and seashell
PT4	E6	1–2	Poorly sorted	Gradational	Relatively sharp	2.5Y 2.5/1	Presence of gravel, angular fragments of seashell, seashell and fish bones

* 40 cm Rip-up clast.

).

The three other sediment deposits E2, E4 and E6 were composed of smaller grains, and according to the Munsell color scale, these fine sediments (medium grain size: 0.039–0.630 mm) were black (2.5Y 2.5/1) (Figure 3).

The unit E1 (30–50 cm) is dominated by very coarse and coarse sand (medium grain size: 0.550–0.800 mm), moderately sorted of dark gray color (2.5Y 4/1), with the presence of marine calcareous remains and stones (>3 mm/gravel). This deposit includes two rip-up clasts at 30 and 40 cm of depth; the upper was made of dark gray clay (2.2Y 2.5/1), while the other consisted of gray clay. These rip-up clasts are fine sediments (clay or/and silt) with very poorly sorted, that are formed due to the erosive current flows and are transported some distance. The unit E3 (20–23 cm depth) was made up of medium and fine sand (medium grain size: 0.260–0.360 mm), sorted and with a low presence of a few fragments of seashells. They are observed at the upper and basal relatively sharp contact line (Figure 3 and Figure 4A; Table 3). Finally, unit E5 was formed by gray (2.5Y 4/1) coarse sand and medium sand (medium grain size: 0.400–0.500 mm), these are moderately well sorted. They are observed at the upper relatively sharp and basal sharp contact line with the presence of marine calcareous remains and stones (20 mm/gravel) (Table 3, Figure 3 and Figure 4A).

The E2, E4 and E6 deposits were generally characterized by a medium skewness (−0.8 to 2; Fine Skewed to Very fine skewed) and kurtosis (3 to 7; Leptokurtic), rich in organic matter, plant remains, roots and wood charcoal (Figure 3 and Figure 4A). The deposit E2 (21–23 cm depth) was characterized by finer sediments (0.250–0.330 mm), poorly sorted and a low bulk density (~1.17 g/cm3). They are observed at the upper relatively sharp and basal sharp contact line. The layer E4 (10–11 cm) was characterized by finer sediments (0.140–0.200 mm) and the lowest bulk density (<1.3 g/cm³), with a poorly sorted and upper sharp and basal sharp contact line. This deposit had the highest content of wood charcoal particles. Finally, the E6 (2–3 cm depth) had fine sands (~0.300 mm), moderately well sorted and higher bulk density (1.5–1.6 g/cm³). They are observed at the upper gradational and basal relatively sharp contact line (Figure 3 and Figure 4A; Table 3).

Biogenic lagoonal deposits are not always present. These deposits are characterized by a mean grain of 0.300 mm, a decrease in density, kurtosis from the base upwards (~8; Very leptokurtic) and skewness (>2; Very fine skewed). In addition, they are observed at an erosive contact line between deposits (Figure 3 and Figure 4).

The core PT4 presented the same deposits succession as core PT1, but the deposits were composed of finer grains size than the PT1 core; the range size was 0.180–0.500 mm compared to 0.140–0.800 mm in PT1. The deposits E1 and E5 in PT4 core were characterized by medium and coarse sands (>0.280 mm), poorly sorted (1–2), medium skewness (2–3; Very fine skewed3) and high kurtosis (8 to 22; Very leptokurtic), except in the range of 13–15 cm, in which the values increased skewness (>4; Very fine skewed) and the kurtosis (>28; Very leptokurtic). E2 and E4 deposits were characterized by finer sediments (<0.200 mm), poorly sorted (1.5 to 2.5), high skewness (1.3 to 2.8; Very fine skewed) and kurtosis (4 to 7; Leptokurtic) (Figure 3 and Figure 4B; Table 3).

4.2. Geochemistry

The sediments of Pachingo wetland, are formed by the particles of the surrounding landmasses; namely, the old Cenozoic basin sediment and/or that of the coastal sandy dune. These potential source areas are characterized by the different geochemical compositions used to identify the source of lithogenic particles in the sediment cores. We used Principal Component Analysis (PCA) to trace the origin of sediments. The PCA varimax rotation permitted us to identify two components that explained 80% of the total variance (Figure 5).

Factor 1 accounts for 64.48% of the total variance. It is characterized by high positive loadings for Fe, Ti, Zn and Ni, whereas Ca and Sr have negative loadings. Fe and Ti indicate the dominance of alumino-silicate minerals, and the granulometric distribution related to these elements was in the clay and silt range (Figure 4 and Figure 6). These elements were linked to the fluvial source. The high percentages of Ca and Sr were related to the significant presence of biogenic material and shell debris. The granulometric distribution of sediment samples associated with these elements was the sand fraction (Figure 4 and Figure 6) linked to the marine source. Thus, we use the Fe/Sr ratio to discriminate between the two potential source areas; that is, a high Fe/Sr ratio value indicates a high concentration of fluvial terrigenous particles, and a lower Fe/Sr ratio indicates a higher relative contribution from the coastal sandy dunes. Except in E3, where the geochemical values did not show significant variation.

In the PT1 core, E1, E3 and E5 deposits were characterized by a high content of Sr with values above 600 ppm in E1 and E5. These deposits showed low Fe concentrations (<25,000 ppm), resulting in a low Fe/Sr ratio. An anomaly was found in the E1 unit (40 cm) with a low Sr and high Fe content and, therefore a high Fe/Sr ratio, related to the presence of a rip-up clast (Figure 6A).

Similarly, in core PT4, the deposit E1 presented similar Sr values (~600 ppm), increasing from E5 to E6. Contrarily, the Fe in deposit E2 has higher concentrations than those observed in the fine sediments of the core PT1. The Fe/Sr ratio was highest in E2, with its minimum in E1 (Figure 6B).

The δ¹³C and C/N distributions were used to track the sources of organic matter. The continental plants have a relatively narrow δ¹³C range from −18 to −35‰ [67,68]. In PT1, the highest values of δ¹³C were in the first 4 cm, with values around −20 to −22‰ and C/N values of >70, tending to have a slight increase in units E2 and E4. These values decreased to lower than −26‰ and C/N values of <20 towards the core base. Still, they were lower also in units E1, E3 and E5 (Figure 6).

4.3. Biogenic Components

We identified 32 species taxa of diatoms for PT1 (Appendix A and Appendix B). In PT1, the total contents of diatoms varied between 0 and 1.2 × 10⁶ valves g⁻¹, and the total pollen count was between 0 and 304. We noted that the contents of diatoms and pollen in deposits E2 and E4 were maximum (1.2 × 10⁶, 612; respectively) and well preserved in the finer grain sediments. On the contrary, in deposits E1, E3 and E5, the diatoms and pollen abundances were very low and poorly preserved (<200), with a high content of fragments of diatoms coincident with the coarser-grained deposits (Figure 6; Appendix B and Appendix C).

4.4. Marine and Flood Layer Identification

The combined approach of sedimentology and biogeochemistry analyses helped us to identify six events: the E1, E3 and E5 deposits were clearly associated with marine events according to the Sr and Fe content and the low Fe/Sr ratio. In addition, these deposits had many seashells and stones with a size above 0.200 mm (gravel) and were characterized by low contents of diatoms with high amounts of fragments and low pollen counts. (Figure 3, Figure 4 and Figure 6).

In contrast, three fine deposits (E2, E4 and E6) had a high content of Fe, low content of Sr and thus a high Fe/Sr ratio and more freshwater diatoms (Appendix B). These deposits also presented high content of organic matter and high content of plant remains, roots and wood charcoals. These deposits were associated with fluvial events (Figure 3, Figure 4 and Figure 6).

4.5. Geochronology

²¹⁰Pbxs and ¹³⁷Cs activity profiles of the core PT1 are shown in Figure 7. The maximum ²¹⁰Pbxs activity was observed at 6 cm depth in the core (1.5 dpm/g) and decreased to reach constant activities below ~20 cm depth (~0.46 dpm/g). The ²¹⁰Pb activity was variable with depth, which was caused by differences in the grain size and bulk density; therefore, radioactive decay cannot explain the profile by itself. Thus, the activity within E4, E5 and E6 units is linked to the grain-size variations and sediment reworking [83]. In agreement with this, we assumed that coarse sediments with very low ²¹⁰Pb activities at depths of 3, 5–8, 10–11 and 16–20 cm represented anomalies in these reworked deposits (E4, E5, E6) and therefore were not considered in the age calculations.

For better age estimation, we used ¹³⁷Cs activity. This profile presented a sharp peak (0.08 mBq/g) at 11–12 cm and near to zero above and below this depth (Figure 7). The most common dating method based on ¹³⁷Cs data [84] assumes that depths with ¹³⁷Cs activities in the sediment correspond to the period of atmospheric production due to nuclear bomb testing in 1950–1965; therefore, we considered that the sharp peak at 11–12 cm in the core PT1 was associated with the 1950–1965 period. This result is in good agreement with ²¹⁰Pbxs (1976 ± 7) and confirms our assumptions and the estimated age–depth relationship.

For estuarine shells, the ages were calibrated considering the three known reservoir ages reported for the area, 165 ± 107 years [85], 625 ± 46 years [86] and 442 ± 2 years [87]. We obtained an age range because we could not constrain these different reservoir ages. Thus, we obtained a maximum and minimum age for our ¹⁴C dates (Table 1). The ages obtained at other depths (8, 23 and 36 cm) presented inversions in reworked deposits, complicating the age–depth modeling; therefore, we only consider the ages from wood charcoal data at 23 cm and the shell data at 49 cm to date the E1 and E2 units, which resulted in an age between AD 1468 and 858 (minimal age) at 49 cm and 1664 ± 28 at 23 cm (Figure 7 and Table 1).

4.6. SWAN Model

The results of the ocean wave propagation are shown in Figure 8, with the significant height and direction fields for the 100 × 100 (m) (mid-resolution) and 20 × 20 (m) (high-resolution) domains, respectively. For the 5 July 1984 storm, waves coming from the third quadrant (SW) were caused by the presence of Lengua de Vaca Point, which generated a diffraction pole to the incident wave (Figure 8A). This induced attenuation of the wave height within Tongoy Bay, in which the Hs values were between 1.0 and 2.0 m (Figure 8B).

On the contrary, the event of 10 July affected the fourth quadrant (NW), the sector towards which the Bay of Tongoy is exposed (Figure 8C). This means that the waves did not undergo major transformation processes in their propagation, presenting Hs values above 4.0 m within the bay (Figure 8D).

5. Discussion

5.1. Extreme Events Identification

Core PT1 is located in the Pachingo wetland at ~850 m and 920 m from the coastline—an area sometimes flooded by storms, tsunamis or heavy rain events (Figure 8). The whole core comprises alternating biogenic lagoonal sediment, fine-grained pluvial flood facies or coarse-grained tsunami or storm facies (Figure 3 and Figure 4). The pluvial flood facies are characterized by a high content of organic carbon with a low lithic fragment content (<20%). In core PT1, this facies is usually represented by a grain size that is typically ~3mm (Figure 4). Extreme sea events are characterized by the presence of gravel, fragments of seashells, high density, low skewness and kurtosis (Figure 4), in addition to the high concentration of Sr and almost no presence of diatoms and pollen (Figure 6A). The extreme pluvial flood events are characterized by high concentration Fe, TOC and better conservation and presence of pollen and diatoms (Figure 6B), in addition to a low concentration of Sr and absence of seashells. Biogenic lagoonal facies are not always present (Figure 6). This may be explained by the impact of high-energy events that erode this deposit. In that case, we observe an erosive contact line between deposits. It is also possible that this in-situ biological contribution occurs at a very low sedimentation rate and that it is not visible. In this case, we observe a succession of intense events deposits.

5.2. Reconstitution of Marine Submersion Deposits

The sedimentological, biological and geochemical results described above allow us to identify three submersion events (storms and tsunamis). The three layers identified as marine events E1, E3 and E5 are characterized by an increase in Sr content, by the presence of numerous seashell remains and an increase in grain size (>0.300 mm) except at E3. These three deposits are also characterized by a low content and poor preservation of diatoms and pollen. In E1, the grain size is particularly larger and poorly selected (>0.400 mm). The coarse sand and very coarse sand comprise an important percentage of the sediment composition (0.800 mm, >70%). In this deposit, the highest percentage of diatoms were freshwater or estuarine, these are better preserved marine diatoms (Appendix B). This type of distribution may be due to the fact that the tsunami wave enters with great energy, producing poor preservation of marine diatoms, on the contrary, the return wave has less energy and drags freshwater diatoms with it towards the coast. These high-energy events produce a dilution effect, causing the removal of previous deposits (Figure 4 and Figure 6).

Site sensitivity to overwash deposits may result from different factors such as a change in sea level [39,47,48]. An increase in sea level induces a shift of the barrier or sandy beach landward. Therefore, an increase in sand layers in a sediment core may be due to a sea level change. In our area, the sea level has remained more or less constant during the last 500 years (<0.5 m) [88,89]. Moreover, according to Ota and Paskoff (1993), the position of the sandy beach has not shifted significantly during the last 500 years. They suggest a constant progradation rate of ca. 0.14 m/year, representing 70 m seaward in 500 years. This distance is negligible compared to the core site position of our cores, extracted at 850 m and 920 m from the coastline. Therefore, the minor relative sea-level fluctuations and shoreline changes during the last 500 years did not drastically alter the lagoon’s depositional environment; hence the sensitivity of the site in recording intense paleostorms and paleotsunamis. These three deposits (E5, E3 and E1) are thicker in PT1 and thinner towards PT4; if it is more, E3 disappears. The thickness of these two deposits (E5 and E1), decreases from the coastal area toward the watershed; they are overwash layers—that is, coming from marine incursions during intense sea events (Figure 9). In addition, it has been shown that there is a positive relationship between the intensity of marine submersion and the size of the overwash sand body [39,49,63]. In that case, we may suggest that these three sand layers recorded in the Pachingo’s wetland at more than ~1 km from the coastline were probably formed by three extreme sea events.

The E1 unit has the highest Sr concentrations, which is related to the coarsest grain sizes in the whole core and the significant content of seashells. This deposit also presents stones larger than 3 cm (gravel) and two rip-up clasts at depths of 30 and 40 cm (Figure 4 and Figure 6), suggesting a high-intensity event such as a tsunami [89,90,91,92]. According to the ages estimated in unit E1, this could have occurred between 858 and 1468 AD; in this period, two tsunamis were recorded based on textual and geological archives [93,94,95,96,97,98] Trans-oceanic events occurred in 1361 and 1420. The event of 1361 was generated by an earthquake greater than 8.5 Mw in Japan that affected the Prefectures of Tokushima, Osaka, Wakayama and Nara, and the Awaji Island. Large amounts of damage in the Kii Peninsula of Japan, the destruction of ~1700 homes and more than 60 deaths were caused by the tsunami that occurred in Awa [93,94,95] (NCEI/WDS). According to Satake et al. (2020), this event could have had the required characteristics to affect Chile at that time; however, it has not left strong evidence in sedimentary records in Chile [96]. The second event, which occurred in 1420, was generated by an earthquake of 8.8–9.4 Mw in Atacama (Chile), 354 km further north than our study area. As described by Abad et al. (2020), this may have been one of the most significant events in a supercycle of earthquakes in the Atacama region that caused a rupture of greater than 600 km and probably produced a trans-oceanic tsunami that affected the area of Japan (the 1420 Oei orphan tsunami). Abad et al. (2020) described a large boulder produced by a flood that may have had a height of >18.5 m above high tide level and an inland penetration greater than 284 m from the cliff edge of Cisne Bay in the southern Atacama.

According to this presented background, our E1 deposit should correspond to the Atacama Tsunami (1420), due to its proximity to the place of the epicenter and the magnitude described by Abad et al. (2020).

The following marine deposit, E3, is characterized by a decrease in biological parameters such as TOC, total pollen and Fe concentrations. This deposit is not well constrained by a change in granulometry, which suggests that event E3 was probably of lower energy than event E1 (Figure 4). This deposit is dated 1918 ± 26 AD and could be related to the 1922 tsunami, produced by the earthquake which occurred in Vallenar (central Chile), with a magnitude of 8.5 Mw [20,99]. This event produced a destructive tsunami that affected a large part of the area from Chañaral to Coquimbo. The maximum variation in sea level was recorded at 9 m in Chañaral, 5 and 7 m above high tide in Caldera and in Coquimbo Bay, and the run-up exceeded 7 m in the southern part of the Coquimbo Bay [20,32,100,101,102]. Interestingly, this deposit has similar sedimentological characteristics to the sandy layer associated with the tsunami of 1922 described by DePaolis et al. (2021), which also supports its association with this tsunami event.

The last marine deposit, E5, is characterized by a clear increase in grain size and the presence of seashells, in contrast to the E3 deposit. These characteristics were preserved in both cores PT1 and PT4 cores, indicating a more significant event. This deposit is dated 1986 ± 5 AD and could be associated with the tsunami that impacted Valparaiso in 1985, 312 km further south of our study area. It caused a high tide in Valparaiso with a maximum wave height of 1.15 m. However, this event has been characterized as a low-energy event [103] (SHOA, access date: 15 October 2021 http://www.shoa.cl/s3/servicios/tsunami/data/tsunamis_historico.pdf). Alternatively, the E5 deposit could correspond to the tsunami that impacted Valdivia in 1960, 1085 km further south than our study area. This event hit Coquimbo with a wave height of 2.2 m (SHOA, access date: 15 October 2021 http://www.shoa.cl/s3/servicios/tsunami/data/tsunamis_historico.pdf). Apparently, no significant damages were recorded in our study area due to the orientation of Tongoy Bay. Due to their low energy and direction, we suggest that these two events did not significantly impact our study area. During this period, three high-intensity wave storm surge events affected the Coquimbo region in 1965 and two during 1984; both occurred during an El Niño event. This phenomenon induces a rise in the sea level and produces an intensification of sea events [26,103,104,105]. The first occurred on 10 August 1965 (category V) and had extraordinary features that caused great destruction to coastal infrastructure and the flooding of the coast [24,25,26]. Several regions of Chile were affected, between 21°28′–37°00′ S [24], but Tongoy Bay did not report damages due to the orientation of the bay. The second occurred on 5 July 1984 (category II). This event affected the north and center of Chile. It produced damages to fishing vessels, and the ports were closed as a preventative measure. Both events originated from the SW direction (third quadrant), and Tongoy Bay did not report damages due to Lengua de Vaca Point. This induced the attenuation of the wave height within the bay and the study area (Figure 8A,B). Finally, the third event occurred on 10 July 1984; in contrast to earlier events, it occurred in the NW direction (fourth quadrant). Therefore, the bay was exposed to the waves hitting the coast with Hs values above 4.0 [m] (Figure 8C,D). This event, of category IV, produced boat strandings, port closures and the flooding of the coastal area. This storm caused enormous damages and destruction to the coastal infrastructure of the central coast of Chile, similar to the 2015 event [26]. These storm waves were more intense than the 2015 tsunami [26], and in Coquimbo, the run-up of the storm waves exceeded 3 m. This event was probably strong enough to produce sand deposits more than 1 km from the coastline (unit 5).

5.3. Reconstitution of Flooding Associated Extreme Rain

It should be noted that Pachingo only records extreme flood events. The creek is fed by underwater flow, and only an important superficial flow is seen during heavy rains [62], allowing the transport of sediments from the basin to the lagoon. This occurs in periods of the warm phase of ENSO (El Niño phenomenon) [56], Appendix D.

The three events corresponding with the deposits in units E2, E4 and E6 are characterized by a high content of Fe and TOC concomitantly with a decrease in the grain size, density and Sr content (Figure 4 and Figure 6). These deposits are thicker at PT4 and become thinner at the PT1 core (Figure 9). The E2 deposit has a high Fe content, high Fe/Sr ratio and TOC with a large amount of wood charcoal and roots (Figure 6). These characteristics are also observed in PT4. This deposit is dated 1664 ± 562 AD and could correspond to a humid phase developed between 1550 and 1650 AD, related to the Little Ice Age. This period was characterized by a southern displacement of ITCZ from the mean position, increasing the precipitation southward [106], it is reported in a study based on tree rings reconstruction [107] and the public database SADA (South American Drought Atlas, https://sada.cr2.cl/?fbclid=IwAR0fiTvBJWH76A-pQjQKS3WLyEti68GOZpROKvIb2EnesSGfotp63NzEiHQ, access date: 15 October 2021) [108]. The increase of intense heavy rain and floods during this period strongly affected the rainfall distribution in Peru and Argentina [109,110,111].

The deposit E4 is characterized by a decrease in grain size and an increase in the Fe and the Fe/Sr ratio (Figure 4); it presents a high diatom and pollen content, roots and better preservation of these biological components (Figure 6A). These characteristics are also observed in the two cores. This unit corresponds to 1976 ± 7 AD and is consistent with the flood of 1957 or 1972 related with the El Niño event [56] (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php10; access date: 20 October 2021). (i) According to DGA (General Direction of Waters, Chile) data (ID 04552002) from the Chile Tower station (30.6164° S, 71.3742° W), heavy rains occurred in May 1957 with a value of 229 mm per month (https://explorador.cr2.cl/; access date: 5 October 2021). These intensive rains began on 18–19 May and affected the north and the center of Chile, leaving economic losses that exceeded 8000 million pesos, 20 people dead and more than 4 thousand victims. The Elqui River overflowed, washing away houses and neighboring buildings and breaking electricity and drinking water services [27,28]. This event was classified as the most catastrophic mud flow in the region’s history, and La Serena and Coquimbo were isolated for many days [56]. (ii) According to Ortega et al. (2012), the rains that occurred on 24 August 1972 (occurred during an El Niño event) produced flooding and landslides. There were trees that fell down, power lines were cut, traffic interrupted and the streets of La Serena city seemed like rivers. In addition, the was intense snowfall at Elqui Valley, registered at more than 4 m. According to the official journal, this event was the most intense and persistent since the mud flow of 1957.

Finally, the age of event E6 is characterized by high Fe content, Fe/Sr ratio, C/N ratio and the presence of roots (Figure 4). These characteristics were observed in both cores PT1 and PT4. This unit is dated at 1999 ± 4 AD, and fits reasonably with the 1997 extreme pluvial flooding, concomitant with the strongest El Niño in the last decades. This event caused damages from Puerto Montt (42° S) to Coquimbo (30° S), and the latter was the most affected region, with monthly rainfall that exceeded 150 mm (June and August, with 164.5 and 143 mm, respectively). This event caused significant economic losses due to bridge cuts, electrical lighting cuts and the destruction of houses and buildings near the riverbed and streams, leaving thousands of people affected and isolated [27,28].

Our sediment cores have recorded three flood events between 1500 and 2000 AD, related to El Niño events. Other El Niño and intense flood events occurred between 1500 and 2000 AD, including 1641, 1700, 1840, 1965, 1972, 1983, 1987 and 1991–1992, but these events were not recorded in our sediment cores; most probably, they were eroded by marine events. Otherwise, for the wetland, a flood deposit should not always be associated with a single event, but rather with several heavy precipitation events occurring in a short time. It is a small basin covering 487 km², with a temporary runoff strongly affected during abundant rain forming small lagoons that do not reach the sea. Thus, the Pachingo wetland is not an ideal area to reconstruct the recurrence of flood events [27,28,56,104,106,112]; however, we could recover evidence of the most extreme events that have impacted the region.

6. Conclusions

The cores retrieved from Pachingo wetland provided sedimentological and age information of the last 1000 years. We identify six extreme events related to marine and continental deposits; the deposits named units E1, E3 and E5 were identified as marine events, and E6, E4 and E2 as continental. The E1 and E5 units correspond to tsunamis that occurred before 1500 AD and in 1922, respectively. The E1 unit would be related to the Atacama tsunami recorded in 1420. The E3 unit does not seem to be associated with a tsunami event but with an extreme storm event on 10 July 1984.

Regarding the continental events, the most recent (E6) seems to correspond to the heavy rain event of 1997, intensified by the warm phase of ENSO. The event E4 could correspond to the 1957 heavy rain event, which produced large overflows of rivers. The unit E2, dated at 1664 ± 28 AD, was of greater intensity, producing a significant carry-over of organic material, and could correspond to a humid phase between 1550 AD and 1650 AD.

The six events identified in the sedimentary records agree very well with the main extreme events reported for the zone; unfortunately, the lack of old records (>1930) makes it difficult to establish the events of unit 1 with confidence, but it certainly is a marine extreme event. This information highlights the relevance of these studies that record the sensitivity of these coastal areas to the impact of these extreme events related to the meteorological variability in the SE Pacific. Additionally, it will enrich the Chilean database of past marine and continental extreme events and assist coastal managers in improving coastal zone management and developing adaptation strategies to protect the population from future damage caused by them.