*Article* **Space–Time Variations of the Apparent Resistivity Associated with Seismic Activity by Using 1D-Magnetotelluric (MT) Data in the Central Part of Colombia (South America)**

**Carlos A. Vargas, J. Sebastian Gomez, Juan J. Gomez, Juan M. Solano and Alexander Caneva \***

Department of Geosciences, Universidad Nacional de Colombia at Bogota, Bogota 111321, Colombia **\*** Correspondence: acanevar@unal.edu.co

**Abstract:** In this work, we apply multi-temporal 1D-magnetotelluric (MT) surveys to estimate the space–time variations of the apparent resistivity *ρ<sup>a</sup>* and correlate these changes with seismic activity in the central part of Colombia (South America). We use the time series of the Earth's natural electric and magnetic fields registered at two MT stations of the National University of Colombia Seismological Network (RSUNAL), located in the Eastern Andean Cordillera, in the central part of Colombia, over several days. Assuming that large earthquakes may generate these types of anomalies, we identified positive results for the Mesetas earthquake (Mw6.0, Lon = 74.184◦ W, Lat = 3.462◦ N, H = 13 km-depth, 24 December 2019, UTC 19:03:55), with anomalies registered eight hours before the mainshock. The depth at which the resistivity anomaly was identified coincides with the depth of the earthquake hypocenter. The origin of these anomalies may be associated with the migration of fluids due to the change in the stress regime before, during, and after the earthquake. We hypothesize that before the occurrence of an earthquake, the stress field generates pore pressure gradients, promoting alterations in fluid migration that change the resistivity of the upper crust.

**Keywords:** apparent resistivity; correlation; earthquakes; magnetotellurics; electromagnetic signals; electromagnetic anomalies

#### **1. Introduction**

Variations in the physical properties of the subsoil, such as electrical resistivity, caused by the occurrence of seismic activity have been studied in recent decades. Authors such as Du et al. [1,2] established a non-linear relationship between the magnitude of the earthquakes and the amplitude–duration of the recorded anomaly in the apparent resistivity *ρa*. Additionally, various electromagnetic (EM) anomalies have been reported, both in the ionosphere and subsurface, for large-magnitude seismic events [3,4]. These mentioned studies use the maximum entropy method (MEM) and spectral density to analyze their signals.

This study seeks to estimate the variation, both in time and in space (depth), of the *ρ<sup>a</sup>* of the subsoil using seismic events of magnitude > 4.5 in Colombia, using the magnetotelluric (MT) method. In the MT method, deepened and applied to geophysical exploration by Cagniard [5], relationships are obtained that allow the electrical resistivity and its distribution in the subsoil to be inferred from the relationship between the fluctuations in the orthogonal components of the natural electric and magnetic fields. These equations were coupled in an algorithm developed in Python to obtain *ρ<sup>a</sup>* variations as a function of time from the time series of the registered electric and magnetic fields. The method allows us to reach depths from tens of meters to hundreds of kilometers [6], considering that this depends on the used time window. By choosing different time windows, we can achieve a depth ranging between ~100 m and ~100 km. It is possible to reach greater depths, for which it would be sufficient to vary the time window used.

**Citation:** Vargas, C.A.; Gomez, J.S.; Gomez, J.J.; Solano, J.M.; Caneva, A. Space–Time Variations of the Apparent Resistivity Associated with Seismic Activity by Using 1D-Magnetotelluric (MT) Data in the Central Part of Colombia (South America). *Appl. Sci.* **2023**, *13*, 1737. https://doi.org/10.3390/ app13031737

Academic Editors: Eleftheria E. Papadimitriou and Alexey Zavyalov

Received: 20 December 2022 Revised: 19 January 2023 Accepted: 27 January 2023 Published: 29 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. Materials and Methods**

#### *2.1. Electromagnetic Signals of Seismic Origin*

The scientific community has studied the possible relationship between electromagnetic phenomena and seismic events since the second half of the 20th century. There is more and more evidence about a connection between the generation of electromagnetic signals and the occurrence of seismic events (in addition to the well-known variations of the Earth's magnetic field associated with the incidence of the solar wind). This field of study requires combining various disciplines, such as physics, geology, and even biology [7,8]. Anomalies in electromagnetic signals associated with strong seismic events have been identified, and physical mechanisms have been proposed to explain them [9–13].

The study of seismic–electromagnetic phenomena has acquired greater importance since the 1980s [7]. Anomalies were found in the behavior of electromagnetic signals before some strong earthquakes, shedding new light on the study of seismic precursors. Several authors have proposed different mechanisms for generating anomalous electromagnetic signals associated with seismic activity.

For example, (1) a mechanism that could explain the anomalies is the displacement of the limits between cortical blocks of high and low conductivity [14]. Due to this, the electrical currents in the subsoil would be altered, also altering the induced magnetic fields. (2) Another mechanism could be the electrokinetic effect, which consists of the fact that, in the presence of fluids in the subsoil, the solid phase can be electrically charged, while the liquid phase acquires the opposite charge [13]. The latter, when flowing, generates an electric current that can induce magnetic fields. (3) On the other hand, the piezomagnetic effect occurs when the magnetization of ferromagnetic rocks changes when subjected to stress, giving rise to an alteration in the magnetic field in the form of very low-frequency signals [15]. (4) An alternate mechanism could be associated with the formation of microfractures in the rock, in the presence of fluids. As the pressure increases before the seismic event, the rocks dehydrate and lose conductivity. These ruptures can also cause the separation of ionic bonds, generating a variation in the potential difference, which can generate electromagnetic signals [13,15,16].

Zhang and Shen [17] have proposed four stages before large-magnitude earthquakes based on the Wenchuan event, M8.0, on 12 May 2008. The first stage would consist of the accumulation of mechanical stresses. The pressure would close the subsurface microfractures, altering their *ρa*. In the final phase of this stage, electromagnetic signals would be recorded. In the case of the Wenchuan event, this stage lasted for more than two years. The second stage would consist of a blockade of the subsoil. At this point, the microfractures would have closed entirely, and tectonic activity (and therefore seismicity) would have been reduced to a minimum, avoiding the generation of electromagnetic anomalies. The third stage would consist of an "unblocking" and an expansion of the subsoil. The temperature would rise, generating electromagnetic signals in the infrared. Microfractures would re-form, and small seismic events would occur with very lowfrequency electromagnetic emissions. The fourth stage would take place hours or days before the big event and would give rise to anomalies in the electric field of the ionosphere. Although the authors base their model on the observations associated with this specific event, they could be considered for the study of other events.

Solar activity can also be a source of electromagnetic signals, for which it is necessary to distinguish between these and those possibly associated with seismic activity [18]. Likewise, anthropic activity can generate electromagnetic noise, especially in areas where there are electrical power distribution systems [19]. Other problems associated with the identification of very low-frequency anomalies are their low intensity, their similarity with very low-frequency signals that are not of seismic origin, as well as the difficulty in identifying the focus of the earthquake to which the anomalies may be associated, before the earthquake occurs [14].

#### *2.2. Magnetotelluric Sounding (MT)*

The magnetotelluric method studies the penetration and propagation of electromagnetic waves inside the Earth, associated with the action of electrical storms and/or the incidence of the solar wind on the Earth. The method is based on measurements of the natural electric fields on the Earth's surface (by means of 2 perpendicular electric dipoles, *Ex* and *Ey*) and magnetic (using 3 coils or triaxial magnetometers, *Bx*, *By*, and *Bz*).

Among the advantages of the magnetotelluric method is the fact that it is based on using a natural electromagnetic source. The Earth's surface partially reflects the fluctuating electromagnetic fields originating in the ionosphere, and the ionosphere reflects the returning fields again due to its conductive characteristics. This happens repeatedly, so that the fields eventually have a strong vertical component and can be considered as vertically propagating plane waves, characterized by covering a broad spectrum of frequencies. These fields penetrate the ground and induce telluric (electric) currents, which generate secondary magnetic fields. The telluric currents, detected by two pairs of electrodes that make up each pair of dipoles, are oriented in the N–S and E–W directions. The three components of the magnetic fields are measured: the vertical component and two horizontal components, parallel to each one of the electrical components [6,20].

This method provides information about resistivity (conductivity) values for much greater depths than artificial source induction methods. Using long-period signals in the range from 10 to 1000 s, the MT method is relevant to the investigation of the crustal and upper mantle structure [20].

The two electrical and three magnetic components of an MT sounding are recorded continuously during long observation intervals, which can last hours or days. The registered magnetic fields have external contributions from the ionosphere and internal ones related to the distribution of the induced currents. Although an MT sounding can be carried out in a range of infrasound frequencies (ƒ ≈ 10–100 Hz), its main application is the determination of electrical conductivity (resistivity) values at great depths, using very low frequencies (ƒ ≈ 1 Hz) [6,20].

#### *2.3. Magnetotelluric Fundamentals*

The theoretical principles of the method lie in Maxwell's equations, which describe the behavior of electric and magnetic fields, as well as their interaction.

Based on the assumptions of this method [6,21] and the electromagnetic induction phenomena described by Ampere's and Faraday's laws, we define *skin depth* (*δ*). This concept describes the distance, in-depth, which the electromagnetic wave has traveled when its amplitude has decreased by a factor of 1/e, assuming a homogeneous halfspace. This value depends on the resistivity of the medium (*ρ*) and period (*T*) of the electromagnetic wave. Using SI units, the *skin depth* is approximated as:

$$
\delta \approx 500 \sqrt{\rho T} \,\tag{1}
$$

Remembering that period (*T*), frequency (*f*), and angular frequency (*ω*) are related by *T* = 1/ *f* = 2*π*/*ω*, we can use this last equation to analyze the behavior of the resistivity (*ρ*) as a function of the period (*T*), and thus obtain a first approximation of the depth of investigation. This depth (*δ*) is used to approximate the magnetotelluric data.

#### Transfer Function and Apparent Resistivity *ρ<sup>a</sup>*

The transfer function term involves an Earth model describing a linear system with a predictable input and output. Here, the transfer function *C* of Schmucker–Weidelt [6,22,23] was used, which depends on the frequency and can be calculated from the field measurements in their orthogonal directions—in other words, by measuring the electric field in the E–W component (*Ex*-HQE) and the magnetic field in the N–S component (*By*-HFN), or with their equivalents (*Ey*-HQN, *Bx*-HFE), as follows:

$$\mathcal{C}(\omega) = \frac{E\_x}{i\omega B\_y} = -\frac{E\_y}{i\omega B\_x} \tag{2}$$

Then, we can calculate the *ρa*, defined as the resistivity's average of a homogeneous half-space. Like *C*, the *ρ<sup>a</sup>* is also expressed as a function of frequency [6]:

$$\rho\_d(\omega) = \left| \mathbb{C}(\omega) \right|^2 \mu\_0 \omega\_\prime \tag{3}$$

where *μ*<sup>0</sup> is the magnetic permeability of free space.

In this study, we use the apparent resistivity *ρa*. However, it is essential to note that for MT soundings, being *C*-complex, it is also important to obtain the phase parameter, which, like these previous functions, will also depend on the frequency:

$$\phi\_{1-D}(\omega) = \tan^{-1}\left(E\_x(\omega) / B\_y(\omega)\right) \tag{4}$$

Note that these MT parameters are usually plotted as a function of the period (*T*), as seen in the Results and Discussion Section.

#### *2.4. Kp Index*

According to NOAA's Space Weather Prediction Center [24], this index quantifies the impact on the Earth's horizontal magnetic field measured on the surface due to the geomagnetic activity, that is, the emission of charges and high solar radiation. This represents a significant source of noise for time series and data. The Kp index activity is classified as follows (Table 1):

**Table 1.** Kp index activity classification.


A Kp index > 4 (high solar activity) would imply a high disturbance in the terrestrial magnetic field and insufficient data reliability in the chosen time interval. This value was downloaded from the NOAA (National Oceanic and Atmospheric Administration) database on its website [24].

#### *2.5. Mutiparameter Station*

The USME station is a multi-parameter station that is part of the Seismological Network of the National University of Colombia [25]. Figure 1 shows the distribution of the network's stations in the central region of Colombia. The USME station has a broadband triaxial seismometer, triaxial magnetometer, and non-polarizable electrodes.

**Figure 1.** Distribution of the stations of the Geophysical Network of the National University of Colombia in the central region of Colombia.

#### 2.5.1. Sensors

The sensors used in multi-parameter stations are [26]:


Figure 2 shows the instrumental deployment corresponding to the USME station.

**Figure 2.** Schematic representation of the location of the orthogonal dipoles, corresponding to the USME station.

#### 2.5.2. Acquisition Hardware

The digitizer is based on a Raspberry Pi 3 CPU with 1.2 GHz processing and 2 GB of RAM, with an independent video chip. Different peripherals are integrated into it via GPIO or USB. The 24-bit Analog–Digital Converter (ADC) is connected via serial communication protocol, having eight single channels to sample up to 500 sps that can be expanded to 16 channels. To achieve time synchronization with an error of fewer than 10 μs, we use a GPS connected via asynchronous serial communication, with the possibility of an external high-gain antenna to improve the view of the sky. Data transmission is via 3G modem. The power supply is obtained using a 5 V DC source at 2.5 A that is provided using a solar panel. The maximum consumption is around 10 W. The entire system is installed in a hard box.

#### *2.6. Data Processing*

The MT data processing in this study begins by converting a time series of magnetic and electrical signals to their frequency domain (using the Fourier transform (FT)), a stacking process over a 1-h time window (it is desired to see changes in the resistivity behavior each hour), to later operate their respective components, calculate the Schmucker– Weidelt transfer function from Equation (2), and finally obtain their associated values of *ρa*. This procedure was performed by developing our algorithm code in Python, which can be obtained by consulting the authors. The processing was performed following all the mathematical guidelines given in [6].

#### 2.6.1. Time Windows

To estimate the space–time *ρ<sup>a</sup>* variation, particular time intervals should be defined depending on the changes to be seen. Therefore, we define three principal time windows: *total window*, *partial window*, and *calculation window*.


• Calculation window: It is the time interval over which we will calculate the Fourier Transform (FT) and, therefore, the variable used as the period (*T*) in the MT method equations (e.g., 10 s, 100 s, or 1000 s, depending on the depth of penetration to be obtained).

In this study:


2.6.2. Apparent Resistivity Dispersion—Apparent Resistivity Curve

The *ρ<sup>a</sup>* parameter (*ρa*), using Equations (2) and (3), will result in several dispersion points, as observed in Figure 3. Note how the figure is plotted as a function of the period (*T* = 1/ *f*) for ease of analysis (the longer the period, the greater the *skin depth*). After obtaining all the associated *ρa*, a statistical data treatment is carried out for values of frequencies (periods) chosen in such a way that a reliable resistivity curve is obtained.

**Figure 3.** Apparent resistivity *ρ<sup>a</sup>* dispersion as a function of the period (*T*). The gray dots represent the raw data computed from Equation (3). The red dots represent the fitted curve used for the apparent resistivities of the partial time window ('xy' means HQE-HFN).

To choose the period points to be evaluated, we consider the conditions stipulated in [6]:


Points are selected as follows:

$$\begin{array}{l} f\_1 = f\_{\text{max}}, \\ f\_2 = f\_{\text{max}} / \sqrt{1.8}, \\ f\_3 = f\_{\text{max}} / 1.8, \\ \vdots \\ f\_k = f\_{\text{max}} / 1.8 \sqrt{\left(k - 1\right)}, \end{array} \tag{5}$$

where *fmax* = 10 Hz (sample rate = 20 Hz).

After obtaining the desired frequencies, we observe that the dispersion of the data for each point will have a Gaussian behavior, as seen in Figure 4. This allows us to find a single *ρ<sup>a</sup>* value associated with a single frequency/period and the corresponding standard deviation. Figure 3 shows an example of a curve obtained from fitting our dispersion data.

**Figure 4.** Distribution of the first six frequencies chosen to represent the apparent resistivity curve.

Now that we have the *ρ<sup>a</sup>* curve for each hour, an *ρ<sup>a</sup>* map is made in the total window (days), thus allowing precise observation of the space–time changes where the target seismic event occurred.

#### 2.6.3. Apparent-Resistivity *ρ<sup>a</sup>* Map

The *ρ<sup>a</sup>* map (Figure 5), which is obtained from our algorithm code, can be seen as a matrix of values of *N* rows and *M* columns (Equation (6)) where each column represents a 1-h time window, and each row represents each frequency/period point (*N* = 'Number of point periods' and *M* = 'Number of partial windows'). Each column corresponds to the apparent resistivity-fitted curve (Figure 3); the only difference is the color assigned to each resistivity value (Figure 5).

$$\begin{array}{ccccc}\rho\_{a\_{(0,0)}}\rho\_{a\_{(0,1)}} & \cdots & \rho\_{a\_{(0,M-1)}}\\\rho\_{a\_{(1,0)}}\rho\_{a\_{(1,1)}} & \cdots & \rho\_{a\_{(1,M-1)}}\\\vdots\\\rho\_{a\_{(N-1,0)}}\rho\_{a\_{(N-1,1)}} & \cdots & \rho\_{a\_{(N-1,M-1)}}\end{array} \tag{6}$$

**Figure 5.** Anomaly registered before the 6.0 Mw Mesetas earthquake calculating the apparent resistivity from USME HQE-HFN channels. (Format of digitized data on the time axis: year-monthday hour).

As a second option, moving on to the interpretation of the *ρ<sup>a</sup>* map, this can also be seen as a depth profile. Each period in each column means a different and greater depth; let us remember the relation of Equation (1), which tells us that the longer the period, the greater the *skin depth* (keeping in mind that this parameter also depends on the resistivity of the medium).

#### 2.6.4. Anomaly Map

As the *ρ<sup>a</sup>* color map can show some anomalous behaviors at specific depths (space– time variations), sometimes this phenomenon cannot be seen clearly. Therefore, it is necessary to create an anomaly map that represents the same time interval and depth ranges.

The values of the anomaly map (shown in the lower part of Figure 5) are the result of finding how anomalous each value of the *ρ<sup>a</sup>* map is, that is, quantifying how much each one varies according to the average of the total window. First, the average value of the window is found, which will be a column of *N* rows resulting from the average resistivity for each period:

$$\begin{array}{ccccccccc}\rho\_{\mathfrak{a}\_{(0,0)}} & \rho\_{\mathfrak{a}\_{(0,1)}} & \cdots & \rho\_{\mathfrak{a}\_{(0,M-1)}} & \sum\_{i=0}^{M-1} \rho\_{\mathfrak{a}\_{(0,i)}}/M & \overline{\rho\_{\mathfrak{a}\_{0}}}\\\rho\_{\mathfrak{a}\_{(1,0)}} & \rho\_{\mathfrak{a}\_{(1,1)}} & \cdots & \rho\_{\mathfrak{a}\_{(1,M-1)}} & \sum\_{i=0}^{M-1} \rho\_{\mathfrak{a}\_{(1,i)}}/M & \overline{\rho\_{\mathfrak{a}\_{1}}}\\\vdots & \implies & \vdots & & = & \vdots \end{array} \tag{7}$$

Starting from the matrix, which represents the apparent resistivity map (left part), then (⇒) we average each row (middle part), which is equivalent (=) to an average value corresponding to each row (period); we get:

$$\rho\_{a\_{(N-1,0)}} \rho\_{a\_{(N-1,1)}} \cdots \rho\_{a\_{(N-1,M-1)}} \sum\_{i=0}^{M-1} \rho\_{a\_{(N-1,i)}} / M \overline{\rho\_{a\_{N-1}}}.$$

thus, obtaining an average *ρ<sup>a</sup>* value for each period.

Finally, each value of the initial apparent resistivity matrix (*ρa*(*i*,*j*) ) is operated with its corresponding average (*ρai* ) as follows:

$$Anomaly\_{(i,j)} = A\_{(i,j)} = \frac{\rho\_{a\_{(i,j)}} - \overline{\rho\_{a\_i}}}{\rho\_{a\_i}} \cdot 100,\tag{8}$$

to obtain a percentage of the anomaly. This value means how different each *ρ<sup>a</sup>* value (*ρa*(*i*,*j*) ) is from the average (*ρai* ) of the total window, e.g., getting an anomaly value of 300% means that we have a *ρ<sup>a</sup>* value three times greater than the average.

By obtaining all values of the anomalies, we will be able to plot a map with the same dimensions as the *ρ<sup>a</sup>* map. This map is presented in grayscale (Figure 5) and allows us to see more clearly the area where an anomalous resistivity value occurs according to the average value of the total window.

#### **3. Results and Discussion**

Considering several variables such as solar activity, *ρa*, and associated anomalies, the results of this study provide visual information for interpreting in a simpler and more complete way. Figure 5 shows the space–time changes of the *ρ<sup>a</sup>* for the USME station during the time interval when the Mesetas earthquake occurred (Mw6.0, Lon = 74.184◦ W, Lat = 3.462◦ N, H = 13 km-depth, 24 December 2019, UTC 19:03:55) and its associated aftershocks.

Solar activity (Kp index): Average low solar activity (<2). A behavior without significant variation is observed throughout the time window, except for a small peak of 2.5 on 26 December.

Apparent resistivity: The variation of the *ρ<sup>a</sup>* depending on the depth *H* can be classified into three ranges of periods (Table 2):


**Table 2.** Variations of the *ρ<sup>a</sup>* depending on the depth *H*.

*Associated anomalies:* There is a clear anomaly presented 11–16 h before 24 January (Table 3):

**Table 3.** *ρ<sup>a</sup>* anomaly depending on the depth *H*.


Figure 6 shows the *ρ<sup>a</sup>* curves found for a partial window before, during, and after the anomaly found (black: *ρ<sup>a</sup>* curve before the anomaly, red: *ρ<sup>a</sup>* curve during the anomaly, blue: *ρ<sup>a</sup>* curve after the anomaly). A slight increase can be observed in the periods (depths) involved in the red curve, which means that at this specific time interval, there was a slight increase in the *ρ<sup>a</sup>* of the subsurface between the periods shown in Figure 5. Geologically speaking, and because resistivity in rocks largely depends on their saturation fluids, this anomaly can be related to a momentary migration of fluids in the rock due to stresses produced in the mainshock, and then a return to initial equilibrium.

**Figure 6.** Comparison between three apparent resistivity curves, each corresponding to a different hour interval of December 24 (UTC). The red curve represents the time when the anomaly is observed.

The same procedure was done for the TUNJ station (also RSUNAL station) to rule out the idea that the found anomaly simply meant instrumental noise at the station. The results are illustrated in Figure 7. In this case, we see more anomalies, which are interpreted as external noise due to their periodicity. It is observed that there is also a disturbance of lesser magnitude in a depth range (periods) close to the depth of the seismic event, even if the event occurred at a greater distance from the station (~250 km).

**Figure 7.** Anomalies registered at TUNJ station for the same Mesetas events.

#### **4. Conclusions**

The particularities identified in the case of the earthquakes in Mesetas yield suggestive results of an *ρ<sup>a</sup>* anomaly in the subsoil eight hours before the first recorded seismic event (6.0 Mw), a duration 4–5 h from the start of the anomaly and within a range of periods between 1–30 s (depths between 1.5–40 km). It is also shown that this disturbance does not correspond to instrumental noise since it also exists and is registered at the TUNJ station. The origin of these anomalies may be associated with the migration of fluids due to the change in the stress regime before, during, and after the earthquake. We hypothesize that before the occurrence of an earthquake, the stress field generates pore pressure gradients, promoting alterations in fluid migration that change the resistivity of the upper crust.

The anomalies detected and observed at the USME multi-parameter station do not show a direct and causal relationship with seismic events of less than 5.8 Mw or distances greater than 120 km. However, since the notable anomalies in the designed algorithm must exceed 100% (double the average value), it is not ruled out that there are patterns of changes in the *ρ<sup>a</sup>* of the subsoil of lower percentages related to these events. For these cases, we have used a minor anomaly scale.

**Author Contributions:** Conceptualization, C.A.V., J.S.G., J.J.G., J.M.S. and A.C.; Methodology, C.A.V., J.S.G. and J.J.G.; Software, J.S.G. and J.J.G.; Validation, J.S.G. and J.J.G.; Formal analysis, C.A.V., J.S.G. and A.C.; Writing, C.A.V., J.S.G. and A.C.; All authors have read and agreed to the published version of the manuscript.

**Funding:** Colombia Móvil–TIGO (Colombia) contributed with the data transmission, through mobile network, from remote sites to the data processing center, located at the National University of Colombia at Bogota.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Once the paper is approved, the dataset will be openly available.

**Acknowledgments:** The authors want to express their gratitude to the Universidad Nacional Colombia at Bogotá and Universidad Antonio Nariño at Bogotá. Without their support and mutual collaboration, the permanent operation of the multi-parameter network would not be possible. The authors thank the Special Issue Editors (Guest Editors), for their patience, as well as two anonymous reviewers, for their valuable comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
