Hydrometeorological characterization and estimation of landfill leachate generation in the Eastern Amazon/Brazil

Study sites

Matapi River basin is located in the coastal-estuarine sector of Amapá State –Brazil -, it covers Macapá, Santana and Porto Grande counties (Figs. 1A and 1B). Matapi River is tributary of Amazonas River, in its shallowest stretch. It is influenced by semidiurnal mesotides that often influence pollutants’ hydrodynamics and dispersion in the drainage system, and forest environments in flooded areas’ plains (Cunha et al., 2011). There is small human influence on this basin, although several human activities are evident in stretches close to its mouth: urbanization, livestock, agriculture, industrial district (Cunha et al., 2011; Cunha et al., 2021; Silva, Lima & Tavares-Dias, 2016; Neves & Tavares-Dias, 2019; Abreu et al., 2020; Souza et al., 2021; Felix et al., 2021) and ASMM’s expansion. The study site is inserted in Matapi River basin (0°10′19,63″N; 51°8′37,53″W), mainly in the drainage sub-basin (DSB), whose area is 11.08 km². In geographic terms, the head of this DSB is located the ASMM (Fig. 1C), close to the left side of BR-156 road, 13 km from Macapá’s international Airport, and 3.5 km from the main channel in Matapi River (Leal, 2012).

The region’s topography presents limited natural precipitation drainage capacity, which is characteristic of very plain and winding rivers, at relatively low flow speed (Cunha et al., 2011; Cunha & Sternberg, 2018; Felix et al., 2021). On the short-term, it records speed of ≈1 m s⁻¹ through 2/5 of the complete cycle of tides (12.9h). With respect to other periods-of-time, speed is approximately 0.5 m s⁻¹. Because the study site is close to the equator line, the region has two well-defined climatic seasons: the rainy season (Jan–Jul), with mean precipitation Pm = 2,100 mm, and the dry season (Aug–Dec) (Pm = 178 mm) (Galardo et al., 2009; Kuhn et al., 2010).

The study methodology followed the stages presented in the flowchart (Fig. 2):

Generation of medium flow rates

This topic is about describing the flow generation process during normal weather events. That is, in a typical climatic year of average precipitation expected for DSB. Therefore, we describe the morphometric characterization, based on the analysis of climatic variables of precipitation and, eventually, temperature.

DSB was designed based on the Global Digital Elevation Model (GDEM) Aster 30 × 30. Its morphometric featuring was summarized as follows: (a) area (A) and perimeter (P), (b) form parameters, compactness coefficient (Cc), form factor (Ff), elongation ratio (Re), length and drainage index; (c) relief parameters related to terrain topographic variations, mean slope (j), mean height (Hm); and (d) linear measurements (Lb), drainage density (Dd) (Table 1).

Climate variables such as mean temperature (Tx) and relative humidity (Ur) were calculated through regression test based on the data series of Fazendinha stations –Fz (0°3′S, 51°6′36″W, altitude 14.46 m), Macapá - Mp (0°2′6″N, 51°5′20″W, altitude 17 m) and Porto Grande - PG (041′40″N, 51°24′15″W, altitude 84 m) (1967–2018 data) (Instituto Nacional de Meteorologia, 2000) (Fig. 3B). Wv and Rs were found in the World Water and Climate Atlas of the International Water Management Institute (IWMI), at geographic coordinates 0°10′15″N, 51°8′37″W (centroid ASMM).

Hydrometeorological data were analyzed through the Penman Monteith method, which was modified by FAO (Allen et al., 2006) to generate ETo Eq. (1) in free software CROPWAT v8.0, as well as through processed information on the meteorological variables Fig. S1 and Table S1). (1)${\displaystyle ETo=c\left[w\mathrm{\ast }Rn+\left(1-w\right)\mathrm{\ast }\left(f\left(u\right)\mathrm{\ast }\left(ea\left(1-\frac{Ur\text{\%}}{100}\right)\right)\right)\right]}$ where ETo = Potential Daily Evapotranspiration (mm day⁻¹); c = Penman adjustment factor; w = Penman’s weighing factor; Rn = total liquid radiation (mm day⁻¹); ea = water steam pressure at saturation (mbar); Ur = Relative Humidity; $f\left(u\right)=\left(0.25\ast \left(1+\frac{Wv}{100}\right)\right)$ = wind function; Wv = mean wind speed measured 2.0 m from the ground (km dia⁻¹).

The time series used for hydrometeorological monitoring were obtained from the stations of Fz, Mp, PG (INMET - longest and faulty series) and ASMM (shortest and current series) over the time interval between January 1, 1967, until December 31, 2018. However, possible series inconsistencies may cause significant errors. This fact occurs when gaps in the identification or elimination of data arise (Tarazona, 2005).

The graphic method was adopted to plot hydrograms (Fig. S2) that are visually simpler to show the outcomes (World Meteorological Organization (WMO), 2011; Liu, Wang & Liu, 2017; Moreira et al., 2020), as well as to detect and identify likely drastic changes in the historical series and/or in their trends. Accordingly, it was possible to estimates the periods-of-time whose intervals are scarce (Autori dad Nacional del Agua, ANA; Lema, Plaza & Cecilia, 2019) or reliable (Serrano et al., 2012; Shi et al., 2019).

Student’s t-test assessed whether the different precipitation measurements in the selected periods-of-time were significant (p-value <0.05) by determining likely drastic changes. The F-Fisher test assessed whether the selected periods-of-time presented standard deviations significantly similar or different (p-value <0.05) Eq. (2). Whenever outcomes presented significant differences in the series (Table S2), they were adjusted through Eq. (3) (Villón, 2016), in order to correct the information (Table S3) based on the following statistical impositions: ${\displaystyle \left|Tc\right|\le \left|Tt\right|\left(95\text{\%}\right)\Rightarrow \overline{x_1}=\overline{x_2}\left(\text{means are similar}\right)}$ ${\displaystyle \left|Tc\right|>\left|Tt\right|\left(95\text{\%}\right)\Rightarrow \overline{x_1}\ne \overline{x_2}\left(\text{different means or there are changes}\right)}$ ${\displaystyle Fc\le Ft\left(95\text{\%}\right)\Rightarrow S_1\left(x\right)=S_2\left(x\right)\left(\text{standard deviations are statistically similar}\right)}$ (2)${\displaystyle Fc>Ft\left(95\text{\%}\right)\Rightarrow S_1\left(x\right)\ne S_2\left(x\right)\left(\text{standard deviations are statistically different}\right)}$ (3)${\displaystyle x_{\left(t\right)}=\frac{x_t-\overline{x_1}}{S_1}.S_2+\overline{x_2}}$

where x_(t) is the corrected value, x_t is the value to be corrected, $\overline{x_1}$ is the average of data from the doubtful period, S₁ is the standard deviation of data recorded for the doubtful period, $\overline{x_2}$ is the average of data in the reliable period, and S₂ is the standard deviation of data from the reliable period.

Lack of data shortage in estimated precipitation records obtained through the models were completed by extending data of pre-existing series. Similar analysis was carried out by using several methods (linear regression, multiple regression, means of neighborhoods, mean ratios, correction in comparison to neighbor stations) (National Meteorology and Hydrology Service of Peru (SENAMHI), 2013). The linear regression method consists in establishing a linear regression and correlation between a consistent station and another one that lacks information (Varela, 2007). In order to fulfil and extend the stations (Table S4), they were grouped based on their proximity and similarity to altitude and climate, in case the existence of more representative and reliable stations was confirmed. Accordingly, the series were statistically assessed through the Student’s t and F-Fisher tests, with emphasis on the analysis method known as “double mass” (SearcLUy, Hardison & Langbein, 1960). Correlation, in this method, is found by plotting the value of mean accumulated annual precipitation per station on the abscissa axis, and the accumulated annual precipitation value on the ordinate axis (Fig. S3). The selection of the most reliable model was carried out based on the smallest number of series’ discontinuity in order to generate precipitation per station (Varela, 2007; National Meteorology and Hydrology Service of Peru (SENAMHI), 2013) (Figs. S4 and S5).

Values recorded in the ASMM were taken into consideration to feature precipitation given its effective proximity and representativeness in Thiessen method application. Ratios were altitude similarity, coverage and relief. These stations allowed using the hydrological Lutz Scholz model (Lutz Scholz, 1980), which was defined as deterministic (hydric balance) to calculate monthly runoff for an average year. The Lutz Scholz model was chosen because it allows combining factors that produce and influence generated mean runoffs (Qg) (Sarazú Cotrina, 2018).

The main parameters used in the model were precipitation (Q_p), evaporation (Q_ev), storage (Q_ar) and the natural exhaustion function (Eq. 4). The stochastic mathematical model (Markovian process) was used to calculate Q_g (Eq. 5). Subsequently, the statistical analysis and data correction were carried out through Eq. (3). (4)${\displaystyle a=0.00252\ast Ln\left(A\right)+0.026}$ where a is the exhaustion coefficient per day, A is the DSB area (km²).

This equation enables estimating exhaustion, which was considered median, because it presented retention close to 80 mm year⁻¹, as well as mixed vegetation (grass, woods and cultivated terrains) (Aguirre, 2006). (5)${\displaystyle Q_g=b_1+b_2{Q}_{t-1}+b_{3}P{E}_t+Z.E.\sqrt{1-r^2}}$ where Q_g is runoff in month t; b₁, b₂, b₃ are multiple linear regression coefficients, Q_t−1 are runoff in the previous month, PE_t is the effective precipitation in month t, Z is the often distributed random number (0.1) in month t, E is the standard error of the multiple regression, and r is the multiple correlation coefficient.

Besides generating the Q_g of the study site, additional Q_g were also generated based on area of interest (PMS-01, PMS-02, NSD-01, NSD-02, NSD-03, NSD-04, NSD-05) (Fig. 3C). The PMS-01 is the critical area where the ASMM is located. The criterion for choosing these areas was based on the existence and distribution of small river channels that connect to the mainstream of the DSB. DSB can be defined as exoreic–pluvial classified in the fluvial hierarchy category as third order.

Estimated maximum flow

An information correction process was carried out by multiplying the maximum historical precipitation value within 24 h for each station identified by Fc = 1.13 (Factor of reading correlation, recommended by WMO).

Pmx24h records were analyzed with the Hydrognomon software v.4.01, which allows determining up to 24 distributions per station. This process aimed at assessing the (D_max) adjustment through the Kolmogorov–Smirnov test applied to each distribution (Normal, Log-Normal, Exponential, Gamma, Person III, Log-Person III, Gumbel Max, EV2-Max, Gumbel Min, Weibull, GEV Max, GEV Min, Pareto, Normal L-Moments, Exponential L-Moments, EV1-Max L-Moments, EV2-Máx L-Moments, EV1-Min L-Moments, EV3-Min L-Moments, GEV Max L-Moments, GEV Min L-Moments, Pareto L-Moments L-Moments).

The main aims were to choose the best distribution with the greatest significance (D_max) at p < 0.05 (Kozanis et al., 2010; Güçlü, Şişman & Yeleğen, 2018) and to calculate Pmx24h in short return times (RT): 2, 5, 10, 25, 50 (Marazuela et al., 2019), as well as in longer times, such as the 100 one (Masseran & Safari, 2020). RTs were generated in the Hydrognomon software. The calculated values were interpolated through the Kriging method (Ghosh, Pal & Santra, 2020), with the aid of ArcGIS, to originate the isohyets.

The produced isohyets are curves that connect equal-precipitation points (SearcLUy, Hardison & Langbein, 1960; Herrmann & Bucksch, 2014). As an indirect method, it consists in estimating the Pmx24h in the area of the likely maximum precipitation event in the basin (WMO, 2011). The isohyets were generated from the RT of Pmx24h (2, 5, 10, 25, 50 and 100 years) by taking into account all stations (Marazuela et al., 2019).

Pmx24h of DSB was estimated based on Eq. (6), for each RT. For instance, Pmx24h was found based on the weighed mean relation of the area between isohyet curves and isohyets’ values (Chow, Maidment & Mays, 1994). The rational method (Eq. 7) was adopted to estimate the maximum runoffs (Q_max) (Villón, 2016; Colorado Department of Transportation, 2019). (6)${\displaystyle P=\frac{P_1\ast A_1+P_2\ast A_2+\dots +P_n\ast A_n}{A_{\mathrm{total}}}\mathrm{or}P=\sum \frac{Ai\times Pi}{A_{Total}}}$ where P is Pmx24h of the component (mm), A_i isthe area between each isohyet curve (km²), A_total is the total area of the DSB (km²) and P_i is the mean precipitation between the isohyets (mm) (7)${\displaystyle {\mathrm{Q}}_{\mathrm{g}max}=\frac{CIA}{3.6}}$ where Q_gmax is the maximum estimated runoff (m³ s⁻¹), C is the runoff coefficient –which depends on the vegetal coverage (American Society of Civil Engineers, 2008; World Meteorological Organization (WMO), 2011), on slope and soil type (dimensionless), I is the main maximum intensity –for duration equal to the concentration time –and the given return time (mm h⁻¹) (Ministry of Transport and Communications of Peru, 2011; World Meteorological Organization (WMO), 2011); A is DSB area (km²).

Leaching flow estimates

Predictions from leachate runoff in the projected controlled landfill, which took place in two different ways. The Swiss (Eq. 8) and rational (Eq. 9) methods were applied to estimate leachate production; these methods were substantiated by analyzed basic hydrological parameters (Pm and potential evaporation, data resulting from the application of Lutz Scholz mathematical model). After this stage was over, landfill estimates were summed (for each method) in order to determine leachate runoff, since it would represent the quantitative of the three Accumulation Ponds (AP) in ASMM. Terrain drains headed towards these APs. (8)${\displaystyle Q=\frac{1}{t}P.A.K}$ where Q = mean slurry or leachate flow (L s⁻¹); P = mean annual precipitation (mm); A = landfill surface area –ASMM (m²); t = time in seconds over the year (31,536,000 s year⁻¹); K = coefficient based on the degree of waste compactness (poorly compacted landfills with specific weight of 0.4−0.7 t m⁻³) –estimates show that leachate production ranges from 25% to 50% (K = 0.25 –0.50), which concerns mean annual precipitation and corresponds to the landfill area. (9)${\displaystyle Q_m=\frac{P_m-ES-EP}{t}A}$ where Q_m = mean generated slurry flow (m³ month⁻¹), P_m = mean annual precipitation (mm), EP = mean evapotranspiration (mm); A = contribution area - ASMM (m²); t = time in seconds over the year (31,536,000 s year⁻¹); ES = (P_m x C) = surface runoff (mm); C = surface runoff coefficient (dimensionless runoff).

Predictions helped comparatively analyzing the mean runoffs calculated for different interest areas (PMS-01, PMS-02, NSD-01, NSD-02, NSD-03, NSD-04, NSD-05) (Fig. 3C) and the maximum rainy periods that were projected in different RTs.

Furthermore, considering different RTs for ASMM, the maximum daily leachate production was estimated using Pmx24 h precipitation as an independent variable. In turn, Pmx24 h was calculated from Eq. (8). But it is important highlighting that estimates are necessary because runoff values are currently unknown due to multiple changes and adaptations, such as the expansion of project cells.

Generation of medium flow rates

Based on the Cc = 1.60 calculation, DSB is small and classified as Class III (Oval-oblong to rectangular-oblong) (Viramontes Olivas et al., 2018). It is slightly long (Ff = 0.32) (Delgadillo & Moreno, 2012) and presents low relief (Re = 0.83). The basin’s mean height (H_m) is 11.63 m, and the isometric curves (Fig. 4) are estimates based on the deposition and erosion processes (Hurtrez, Sol & Lucazeau, 1999; Asfaw & Workineh, 2019). In this case, C_l = 5.92 km defines the largest mean extension throughout a straight line between the exit and the water parting. The longitudinal extension of the main channel (E_l = 4.54 km) is defined as the longest effective distance naturally crossed by the water flow in the DSB (Table S5). DSB’s Dd = 0.53 is classified as low density (Delgadillo & Moreno, 2012) and this finding points towards highly resistant strata, dense vegetation and low relief. High D_d values are indicative of weak or impermeable rocks under sparse vegetation and hilly relief (Strahler, 1964; Aladejana & Fagbohun, 2019) (Table S5).

Q_g showed significant variation between monthly means for an average year (p < 0.05) (Table 2). Values 0.071 and 0.723 m³s⁻¹, respectively, represent the runoff in October (minimum) and April (maximum). Similarly, such a process took place in the interest points, PMS-01, PMS-02, NSD-01, NSD-02, NSD-03, NSD-04, NSD-05 (Table 2). Estimates for Q_g differences enabled calculating the concentrations of chemical elements deriving from different sources (soil, vegetation, industrial activity—as observed in the landfill) that interact with the medium (Barros, 2013).

Maximum flow estimate

Table S6 presents the best distribution models and adjustment data (D_max) per station (Fig. 3A). These models allowed calculating Pmx24h, which were conditioned by RT (2, 5, 10, 25, 50 and 100 years). The isohyet curves produced (Fig. 5) for each RT were the final products of the subsequent predictive distribution (Villalta, Bravo deGuenni & Sajo-Castelli, 2020).

Figure 5 represents comparative isohyets of precipitation projected based on the Kriging method in Different RTs. Simulated rain occurrence results with closer values headed Northeast-Southeast (RT = 2 years), and it contrast the other periods from Northeast to Southeast. There are variations in the precipitation map ranging from 94.2 mm (RT = 2 years) to 210.2 mm (RT = 100 years) (Table 3). Variation in isohyets’ direction happens due to changes in wind direction and to increase in precipitation (Pearson III distribution, Table S6, Mp station).

Q_gmax recorded 0.30 of runoff coefficient, and delay coefficient reached 0.60; these values were compatible with basins presenting deciduous forests (Chin, 2013; LMNO Engenharia, Pesquisae Software, 2014) and deciduous foliage (Colorado Department of Transportation, 2019). Total precipitation within one hour (P_d) increased from 46.8 to 104.5 mm and it accounted for increasing Q_gmax from 29.5 to 65.9 m³ s⁻¹, for different RTs (Table 3).

Leachate flow rates

Leachate runoff generation based on the Swiss and rational methods allowed observing significant variation between monthly means recorded for one typical average year (p < 0.05). For instance, values based on both the Swiss (0.0003 and 0.0040 m³ s⁻¹) and rational methods (−0.0055 and 0.0072 m³ s⁻¹) well-represent the generated leachate runoff, mainly in months presenting extreme precipitation—March (maximum) and October (minimum)—since it indicates the lowest values in August and December. On the other hand, values got concentrated below the generated runoff in the interest area (PMS-01), which is highlighted in grey (Fig. 6 and Table S7).

We observed that the variations in the maximum leachate flow estimated by the ASMM in 24 h vary between 28.8 m³s⁻¹ (RT = 2 years) and 64.6 m³s⁻¹ (RT = 100 years). These values are concentrated close to the maximum flow estimated for DSB between 97.5% (RT = 2 years) and 98.1% (RT = 100 years) (Table 4). It is important to note that Q_gmax was calculated disregarding the interference of leachate production.