Observations of Turbulent Heat Fluxes Variability in a Semiarid Coastal Lagoon (Gulf of California)

: Despite the critical role latent (LE) and sensible (H) heat play in turbulent processes and heat exchange in the water–air interface, there is a lack of studies of turbulent fluxes over the surface in semiarid regions. We collected continuous measurements of net radiation ( Rn ) , LE, H, and micrometeorological data at a coastal lagoon in the Gulf of California during 2019 with an eddy covariance (EC) system. We analyzed the time series, considering the North American Monsoon System, the pre-monsoon, and post-monsoon season. Results show that Rn (276 ± 118 W m − 2 ) and turbulent fluxes were higher during the monsoon season (July–September) LE (129 ± 18 W m − 2 ), and H (29 ± 9 W m − 2 ). The monthly average of Rn, LE, and H was highest in June (493.9 W m − 2 ), August (142 W m − 2 ), and May (50 W m − 2 ), respectively. Furthermore, during the monsoon season, the (H + LE )/ Rn ratio (0.74) suggests that more than half of the Rn reaching the coastal lagoon is used for the turbulent exchange of LE and H. During the pre-monsoon, LE (r 2 = 0.36) increases with a higher vapor pressure deficit (VPD), while H (r 2 = 0.66) increases with a higher friction velocity ( u* ) during the monsoon season. Quantitative observations are essential for further research. Z.M.S.-M. Z.M.S.-M.; project Z.M.S.-M.; funding Z.M.S.-M.


Introduction
Water bodies, such as lagoons (i.e., inland and coastal), play an essential role in the energy balance and heat exchange to the atmosphere [1,2]. The exchange of latent (LE) and sensible (H) heat in the atmospheric boundary layer (ABL) along with momentum influences turbulence processes and the development of the ABL [3]. Furthermore, LE links the energy budget and water cycle, providing half of the atmosphere's energy via water vapor [1]. The study of these turbulent fluxes above water surfaces is not new to the scientific community [4][5][6]; however, few studies have aimed to understand these processes in coastal lagoons, especially in a semiarid region.
Turbulence, the flow of mean wind and its perturbations, consists of swirls of motions called eddies, which come in various sizes superimposed on each other [7]. At the surface, the LE flux is the supply that drives atmospheric motion, while the H flux heats the air in the ABL [3]. Over the ocean and coastal surfaces, both fluxes tend to increase with a rise in wind speed [7,8]. In coastal ecosystems, more common than not, LE flux is observed as heat loss by evaporation than gain by condensation [9] and, the H flux is the result of a temperature gradient between the water and the air [10]. Both LE and H move upward from the water surface to the atmosphere [11], that is, the atmosphere gains energy from the water surface, which intensifies atmospheric motions [12].
In general, the vapor pressure gradient will drive the LE flux, and the temperature gradient between the water and the air will drive the H flux [4,13,14]. Furthermore, micrometeorological conditions, such as temperature (T), humidity or the lack of it (vapor pressure deficit, VPD), wind speed (WS) and friction velocity (u*), and solar radiation (SW), as well as, their variability over time (i.e., diurnal cycle and seasons) strongly influence LE and H fluxes [6,[15][16][17]. Most importantly, the amount of solar radiation is vital to both LE and H fluxes, especially during a season with day-to-day fluctuations in cloud cover [5]. Solar radiation is not necessarily a limitation in mid-latitude semiarid regions [18], air temperature (Ta) is high, annual precipitation is low, and potential evaporation is high [19]. In these regions, LE flux over moist surfaces tends to be higher than H flux, in similar proportion over different seasons [20][21][22].
In the semiarid region of northwestern México, the North American Monsoon System (NAMS) a large-scale circulation and rainfall pattern [23], takes place during the summer months from July to September [24]. Although no monsoon season is alike, this synoptic process starts in late June when the Earth's surface warms under intense solar heating [25], hot air rises, and there is a shift in wind patterns, towards the continent (from the South, Southwest) [23]. The result is cloud formation and precipitation in a three month period, which is more than half of the annual precipitation [25].
The turbulent fluxes of LE and H over coastal systems can be studied with different techniques (1) bulk aerodynamic methods upon water-air differences, (2) aerodynamic profile methods, (3) cross-correlation using wind components and air temperature [9], and (4) eddy covariance (EC) [4,22,26]. Worldwide, a broad group of scientists implemented the EC method to measure and monitor energy, water, and carbon fluxes in terrestrial ecosystems [27][28][29]. These measurements are representative of local-scale processes, with the advantage of high-frequency continuous records [30,31]. For continental lakes, marine, and coastal studies, most of the implementation has been in high latitudes [4,[32][33][34][35], with fewer work in semiarid regions [22,36].
The objective of this work is to analyze the variability in a one-year-long dataset, to better understand the net radiation (Rn), latent (LE), and sensible (H) heat of a semiarid coastal lagoon influenced by the North American monsoon. We provide quantitative estimates of (1) the change in the magnitude of Rn, LE, and H, and (2) linkages between turbulent fluxes and micrometeorological drivers. This analysis is critical and advances our current knowledge of water-atmosphere interactions from coastal lagoons, especially in semiarid regions. The results presented here may provide relevant information for future scientific studies.

Site Description
Measurements for this study were performed at the coastal lagoon Estero El Soldado (EES, hereafter) (27°57.248′ N, 110°58.350′ W) from January to December 2019. The EES is a natural protected area (since May 2006) [37] and a Ramsar Site (No. 1982, since 2011) (Figure 1b). The geographical location of this coastal lagoon is in northwestern México, in the central area of the Gulf of California (Figure 1a). This coastal lagoon connects with the ocean via a 40-m-wide and 3-m-deep channel (Figure 1b) [38]. The climate classification is hot-dry defined as BWh [39], based on an adaptation from the Köppen classification system.
The NAMS influences precipitation patterns in this region [40]. Most of the rain falls from July to October, annual precipitation accounts for ~260 mm, while the average annual temperature is ~23 °C with cold winters and hot summers. We obtained these annual estimates from long-term observations (1990-2015) of meteorological station No. 26204 (Empalme,Comisión Nacional del Agua (CONAGUA)-Servicio Meteorológico Nacional (SMN)).
Because the NAMS drives the climatic system in this region, to better understand the Rn, LE, and H fluxes, we considered three seasons in this study: the pre-monsoon, as a transition season from cold winter months to hot summer ones (January to June); the monsoon, when most of the precipitation occurs (July to September); and the post-monsoon, when solar radiation decreases (October to December).

Micrometeorological Measurements
Instruments for eddy covariance (EC) and micrometeorological measurements were attached to a tripod at 2 m above the water (Figure 1d). The tripod was screwed to a floating wood deck to secure it. The deck had a sturdy frame, so it moved up and down with the tide. However, the deck did not move on the horizontal plane; we verified this with an accelerometer (HOBO Pendant G UA-004-64 Data Logger, Onset Computer Corporation, Bourne, MA, USA).
Data collection took place from January to December 2019. Gaps in data were mostly due to energy failure, or maintenance and calibration.
The climatological footprint of the EC system was analyzed using 30-min data for 2019 [48]. Data points coming from 1 to 27 m around the EC system were the most frequent ( Figure 1c). Data beyond 43 m in the horizontal plane were filtered to avoid possible signal contamination from the surrounding vegetation.
Data quality control followed filtering based on [49]; we only used data labeled as "good-0". A second filter was used to avoid outliers beyond ±2σ (standard deviation). Finally, data was lost due to instrument malfunction, maintenance, or power failure. Subsequently, after applying each corresponding filter and accounting for data loss, the total data loss was 29% of the collected data.

Slow Response Measurements
A net radiometer (SN-500, Apogee Instruments Inc., Logan, UT, USA) collected upward (outgoing) and downward (incoming), short (SW) and longwave (LW) radiation at 1.5 m above the sea surface and oriented to the southwest (200°) to avoid solar panel shading. Relative humidity and air temperature were measured using an HMP45 (Vaisala, Helsinki, FIN), and a tipping bucket (TE525, Campbell Scientific, Inc., Logan, Utah, USA) measured precipitation. Data collection was at 1-min frequency. Instruments were connected to a CR1000X (Campbell Scientific, Inc., Logan, Utah, USA) and connected via ethernet to the Smartflux3 through a brainbox (multiple ethernet connector) to integrate 30-min averages.

Tide Data
Tide data was downloaded using the software MAR V1.0 (CICESE, Baja California, México). Daily values include both low and high tide from the Guaymas, Sonora, station; the average of daily values provides information about the spring and neap tides.

Methods
The solar zenith angle (z, °) was calculated as where δ is the declination angle of the Sun, Φ is the latitude, and ω is the hour angle. Radiation and turbulent fluxes above the sea surface were calculated as follow: where radiation components are measured in W m −2 , Rn is net radiation, SW is shortwave, LW is longwave, in refers to downward (incoming) and out to upward (outgoing). The ratio of outgoing SW radiation to incoming SW radiation is referred to as albedo (α) and is calculated as: Latent (LE, W m −2 ) and sensible (H, W m −2 ) heat fluxes are calculated as: where ρ is the density of air (kg m −3 ), cp is the specific heat of air (J kg −1 K −1 ), and Lv is the latent heat of vaporization (J kg −1 ). Furthermore, w is the vertical wind component (m s −1 ), T the air temperature (K), and q is the specific humidity (g kg −1 ), and the primes denote fluctuations with respect to the mean values.
The Bowen ratio (β), is the ratio of sensible heat to latent heat, calculated as: Friction velocity u* (m s −1 ) was calculated as: * = + where u' is the fluctuation in the longitudinal component of the wind velocity, w' is the fluctuation in the vertical wind component, and v' is the fluctuation in the lateral component of the wind velocity.
Vapor pressure deficit (VPD) (kPa) was estimated as: where es is the vapor pressure (kPa) at saturation, and e is the ambient vapor pressure (kPa). Water temperature (Tw) was calculated from the outgoing longwave radiation Stefan-Boltzmann equation where is the outgoing longwave radiation of the water surface, ε is the emissivity of the water surface (0.98), σ is the Stefan-Boltzmann constant and is the water temperature ( ). Therefore, the water temperature can be written as:

Statistical Analysis
The software Matlab® 2015a was used to perform descriptive statistics and probability density function (PDF) analysis. The wind_rose.m file (fileexchange/17748-wind_rose) from the Matlab® file exchange platform was used to analyze wind speed and direction.
A one-way ANOVA with Fisher's least significant difference test was performed to compare the average magnitude of net radiation and turbulent fluxes month to month using Matlab® 2015a.
To analyze the influence of physical factors such as u* and VPD on LE and H turbulent fluxes, data from 10 to 15 h. was used, followed by an ordinary linear regression (OLR) [50,51] using Matlab® 2015a.

Micrometeorological Conditions
At this site in the Gulf of California, the patterns of wind direction and speed vary by season. During the transition season, defined here as the pre-monsoon, the direction from the continent to the ocean (N, NE) is more frequent from January to March; then, from April to June, the wind coming from the ocean (S, SW) increases in frequency (Figure 2a). This pattern of winds from the S and SW remains through all the monsoon season ( Figure 2b). Then during the post-monsoon season (Figure 2c), there is an increase in continental winds from the N and NE. The observed wind patterns are typical of the region, were synoptic process such as the NAMS greatly influence climatic conditions. In Figure 3, we observe the daily accumulated precipitation (mm), and daily averages of wind speed (WS, m s −1 ), air temperature (Ta, °C), vapor pressure deficit (VPD, kPa) and tide height (m), as well as their variability from January to December 2019. Generally, Ta and VPD are lower during the first three months of the pre-monsoon season and the post-monsoon season, and higher during the monsoon season from July to October a typical pattern of this latitude (27° N).
During 2019, annual precipitation reached 258.9 mm, which is close to long-term annual observations for the region. Most of the observed precipitation events recordings were in the monsoon season (152.3 mm). However, some precipitation events occurred during the pre-monsoon season in January (7.62 mm), and March (66.4 mm), due to cold fronts. They also contributed to rainfall events in November (30.2 mm) (Figure 3a).
Based on the annual (2019) measurements, the average wind speed is 2.33 m s −1 (Figure 3b); the range goes from 0.5 to 23 m s −1 when maximum WS is considered (not shown in this figure). However, observations above 10 m s −1 are not frequent at this site. Daily wind speed average during the pre-monsoon and the post-monsoon season was below the annual average and the during the monsoon season above it (Figure 3b).
We observed significant changes in air temperature (Ta) throughout the year; the range goes from 10 to 45 °C (Figure 3c). During the pre-monsoon Ta fluctuates from 10 to 37 °C; during the monsoon season it goes from 27 to 45 °C; then, during the post-monsoon season, it goes from 15 to 33 °C. During the pre-monsoon season, Ta fluctuates the most. The record of mean annual Ta was 25.9 °C. From July (Julian day 200) to September (Julian day 270), Ta was consistently above average, even during the night (Figure 3c). While during the other months, the Ta decreased below average. VPD is highly variable as well (Figure 3d), and ranges from 0.055 kPa to 3.6 kPa; during the premonsoon season, we observe the lowest VPD on average in January (0.6 kPa) and again in December (0.8 kPa) during the post-monsoon, while, during the monsoon season (July-September), VPD increases to about the double especially in August (1.2 kPa). VPD annual mean is 1.023 kPa pre-monsoon, and post-monsoon VPD is below average, whereas during the monsoon, VPD is above average.
The tide in the region is characterized by a semi-diurnal cycle, adding to the dynamic regime of neap and spring tides. Furthermore, during the monsoon months, the water level in the coastal lagoon is above average, and the coastal system is always flooded (Figure 3e). During the pre-monsoon and post-monsoon, the water level is lower-in fact, in some sections of the coastal lagoon, the benthic zone will be exposed to the atmosphere. The channel in the inlet (Figure 1a) will have water all year round, and the minimum height is 0.4 m.

Monthly and Seasonal Variation of Meteorological Conditions, Turbulent Fluxes, and Energy Partitioning
Monthly average wind speed does not vary considerably throughout the year. WS was higher during August (2.6 m s −1 ), July (2.5 m s −1 ), and September (2.4 m s −1 ). The lowest monthly average values occurred during March (2 m s −1 ) and May (2 m s −1 ). The monthly mean friction velocity (u*) variation was similar to WS. We observed the lowest values in January (0.18 m s −1 ) and December (0.22 m s −1 ) during the pre-monsoon and Post-monsoon season, respectively. Furthermore, higher u* records in April (0.27 m s −1 ), May (0.27 m s −1 ), and July (0.28 m s −1 ) ( Table 1).
As mentioned in the section above, VPD and Ta are highest during the monsoon season ( Figure 3). The highest value recorded of VPD was in August (1267 Pa, Table 1); this was also the case for Ta (33.9 °C) and water temperature Tw (32.8 °C) ( Table 1).
Estimates of average Tw are higher than Ta during the pre-monsoon months from January to June, and the post-monsoon months as well ( Table 1). The reason for this is the diurnal variability; Tw does not vary as much as Ta during the night and day because of the heat storage capacity of water. Such diurnal variability, especially from Ta, is less prominent during the monsoon and warmest season of the year (Figure 3c). Table 1. Estero El Soldado coastal lagoon monthly average (and standard deviation in parentheses) of micrometeorological drivers and energy partitioning ratios. From left to right WS (wind speed), u* (friction velocity), VPD (vapor pressure deficit), Ta (air temperature), Tw (water temperature), H (sensible heat flux), LE (latent heat flux), and Rn (net radiation). The SWin was significantly (p < 0.05) lower in December (69.6), and January (87.4), (Figure 4a). SWin was highest in July (481.3), June (350.7), and May (349.8). During the pre-monsoon, SWin was 257.8 ± 107.7 W m −2 , it was highest during the monsoon 313.11 ± 150.9 W m −2 , and lowest during the post-monsoon 132.7 ± 62.7 W m −2 . We observed a significant (p < 0.05) decrease from July to August, associated with cloud formation during the NAMS.

Micrometeorological drivers
Net radiation (Rn) varies through the year, and we observe statistical differences (p < 0.05) at a monthly time scale. Rn was highest in June Latent heat (LE) follows a similar trend to Rn, with high values in August (142.5 W m −2 ), May (138.9 W m −2 ), and April (135.9 W m −2 ). While, significantly lower values (p < 0.05) are in December (73.5 W m −2 ), and January (74.9 W m −2 ) (Figure 4c). LE is highest during the monsoon season 129.8 ± 17.9 W m −2 , followed by the pre-monsoon with 110.6 ± 2.6 W m −2 , and lastly, the post-monsoon season 94 ± 27.25 W m −2 . There is up to a 30% difference between the monsoon and post-monsoon season.
Through this year-round study, LE and H are always positive fluxes, transporting energy from the water surface to the atmosphere. The amount of Rn influences the magnitude of these turbulent fluxes in the system. Furthermore, the energy partitioning through the albedo (α), Bowen (β), and (H + LE)/Rn ratios are a good insight on how much absorbed radiation there is, and how energy is used. Albedo from ocean surfaces without ice is low ~0.08 [6]; therefore, most of the incoming solar radiation is absorbed. In this semiarid coastal lagoon, the seasonal average of albedo was lower in the monsoon (0.07), with the lowest monthly mean value in July (0.03) ( Table 1). During January, February, November, and December ( Figure 3e and Table 1), months from the pre-monsoon and post-monsoon seasons, albedo increases by about ~9%. At the same zenith angle (20 to 40°) from 10 to 15 h, albedo is lower (0.09) with a higher water column (> 0.45 m) than when the water column is lower (< 0.45 m) as, in this case, the average albedo increases (0.11). A shorter water column influences albedo via interactions with the clear sandy sediment and water transparency [52]. Furthermore, ocean surface albedo can be sensitive to zenith angle and wind speed [53], at this site during the same zenith angle (20° to 40°), the albedo was 0.05 either at low WS conditions (< 2.03 m s −1 ) or high wind conditions (> 9 m s −1 ).
Finally, during the monsoon season, we observe that albedo during the same zenith angle (20 to 40°) is higher through cloudy days (0.11) than clear-sky days (0.07). Further research should explore these interactions at different zenith angles and wind speeds.
As expected, the Bowen ratio (β) is highest during the pre-monsoon (0.24) and lowest during the post-monsoon (0.16) season (Table 1). The reduced magnitude of β is the result of the lower values of H. The lowest monthly mean β value was during the post-monsoon season in November (0.08), and the highest was during the pre-monsoon season in May (0.37) ( Table 1).
On the other hand, the ratio of turbulent fluxes to net radiation (H + LE)/Rn also changed with seasons. During the monsoon, the highest seasonal ratio (0.74) was observed, with the highest monthly mean in September (0.80) ( Table 1). The lowest observations recorded were during the post-monsoon (0.54), with December and October being the lowest (0.40). Although the pre-monsoon season average was not the highest (0.58), during January, (H + LE)/Rn ratio value was 0.90. In terms of energy exchange, the annual average of (H + LE)/Rn ratio (0.59) suggests that about half of the Rn reaching the EES coastal lagoon through the year transforms in a turbulent exchange of H and LE. Seasonally, pre-monsoon (H + LE)/Rn ratio suggests that the surface uses only a little more than half of Rn. In contrast, more than half the amount of Rn during the monsoon season (0.74) goes to turbulent fluxes. The fact that there is a higher ratio during the pre-monsoon than during the post-monsoon season relates to energy stored in the water due to heat capacity, that further transfers between season [36,54].
In context with other sites, Rn at higher latitudes is lower than at this coastal lagoon in the Gulf of California (Table 1). In comparison, Rn is about 48% more, LE flux is 40% more, and H flux is 50% more at any given time of the year. With some exceptions, like in Lake W (Canada) during the monsoon, where LE was 44 W m −2 higher than at EES. The higher amount of LE at Lake W links to low solar radiation and a negative flux of H in the bottom of the lake [32]. This phenomenon does not occur at the shallow, and tide-influenced EES. The other exception was the H flux in Ross Barnett Reservoir (USA) [16]; during the post-monsoon, H was about double than in EES (Table 2). Table 2. Net radiation (Rn), latent heat (LE), and sensible heat (H) from water bodies at different latitudes. Data were organized based on seasons in this study, were pre-monsoon (January-June), monsoon (July-September), and post-monsoon (October-December).

Season
Rn

Diurnal Variations of Rn and Turbulent Fluxes
In general, all three fluxes Rn, LE, H, increase at 6 h, reaching a maximum at noon, and decreasing around 19 h ( Figure 5). During the pre-monsoon season, the average values of Rn are similar of those from the monsoon season; however, the variability is much higher during the pre-monsoon season (Figure 5a). Rn during the monsoon drops around noon (Figure 5b), associated with convective clouds, which are more frequent during the NAMS.
From 10 to 13 h, during the pre-monsoon season, there is up to 20% more Rn than in the monsoon and post-monsoon. While at this same time period, LE is 20% more during the monsoon than the pre-monsoon and up to 40% more than in the post-monsoon season. Finally, during the post-monsoon season there is 38% less H flux available, than in the pre-monsoon and monsoon seasons.
It is clear that, with more Rn available in this coastal system, there is an increase in LE and H. Studies focused on better understanding turbulent fluxes over water have reported that the LE always prevails over H, whether is a semiarid [21,22], or a temperate region [55,56]. In this study, we observe how LE consumes almost 40 % of Rn (Table 2 and Figure 5), for all three seasons. The fact that this coastal lagoon is located in the mid-region of the Gulf of California, and therefore influenced by the North American Monsoon System [40], is especially noticeable when analyzing the Rn and SWin.

Micrometeorological Drivers and Turbulent Fluxes
We analyzed the correlation during the pre-monsoon, monsoon, and post-monsoon of vapor pressure deficit (VPD) and latent heat (LE); and friction velocity (u*) and sensible heat (H), using data from 10 to 15 h (zenith angle 20° to 40°).
During the pre-monsoon season, as VPD increases, the turbulent flux of LE increases (r 2 = 0.36) from 50 to 250 W m −2 , in particular from 0.5 kPa to 1.5 kPa at Ta of around 25-30 °C (Figure 6a). In this same season (Figure 6d) as u* ranges from 0.3 to 0.45 m s −1 , H flux increases (r 2 = 0.60) from 20 to 80 W m −2 with Ta between 25-30 °C.
During the monsoon, as VPD increases from 1 kPa to 1.8 kPa, the LE flux increases as well (r 2 = 0.14) from 100 to 300 W m −2 at Ta around 35-40 °C (Figure 6b). Furthermore, when u* increases from 0.2 to 0.6 m s −1 , H flux increases (r 2 = 0.66) from 20 to 120 W m −2 at Ta between 35-40 °C, this is the best correlation for H flux from each season analyzed here (Figure 6e).
Through the post-monsoon season, VPD ranges between 0.5 kPa and 1.5 kPa similar to the premonsoon season, and the LE flux (r 2 = 0.26) goes from 60 to 120 W m −2 with a Ta around 25 °C (Figure 6c). Similarly, to previous seasons, as u* increases from 0.3 to 0.45 m s −1 , there is an increase in the H flux (r 2 = 0.14) of −20 to 80 at Ta around 25 °C (Figure 6f), this correlation, however, is not very strong.
Turbulent fluxes are sensitive to air temperature (Ta), vapor pressure deficit (VPD), and friction velocity (u*), as previously observed in other studies [4,17,22,57]. Interestingly, while H increases with higher u* and high Ta during the monsoon season, LE does not increase as much with higher atmospheric demand during this same period. During the pre-monsoon season, however, VPD can better explain the change in LE flux, especially at higher Ta ( Figure 6). Overall, the sensitivity of turbulent fluxes to micrometeorological conditions changes with the season at this semiarid coastal lagoon. While H flux is more sensitive during the pre-monsoon and monsoon, LE is more sensitive during the pre-monsoon.

Summary and Conclusions
The present study provides observations of Rn, LE, and H flux and its variability at diurnal, monthly, and season scales for an annual period (2019), from a coastal lagoon. The aim of this study was to better understand the turbulent fluxes between the surface water and the air. A limitation of this study, however, is the lack of information on the energy stored in the water column. Despite this limitation, we believe reporting on the Rn and turbulent fluxes of LE and H is important for the atmospheric science, hidrometerology, ecohydrology, and the modeling community. This study is the first to assess Rn, and turbulent fluxes of LE and H in a coastal lagoon of the Gulf of California.
At this coastal lagoon the wind direction after the monsoon season is from the continent to the ocean (N, NE), and it changes towards the continent from the ocean (S, SW) from Abril and through the monsoon season. This air provides moist for cloud formation, SWin and Rn decreased especially in August due to this phenomenon. However, Rn, LE, and H were higher during the monsoon season than any other season. Albedo at this site changes with water column height but not with wind speed at same zenith angle, further research is needed to better understand this dynamic.
The Rn at this site is 40% higher than at other mid-latitudes. Furthermore, LE correlates to VPD during the pre-monsoon and H correlates with u* during the monsoon season. Funding: This research was funded by Mexico's National Council for Science and Technology (Consejo Nacional de Ciencia y Tecnología, CONACYT) with project number CONACYT-278608 and CONACYT-286494, additional funding came from PROFAPI_2018_0033, PROFAPI_2019_0194, and PROFAPI_2020_0020.