A detailed geomorphological map of the studied wetlands was made through stereoscopic photointerpretation of contacts at 1:33000 scale of the aerial photographs of USAF-B flight (1956–1957). The photointerpretation was digitalized and georeferenced using the geographic information system ArcGIS^© (v.10.8.2). The topography of the wetlands basins was analysed by using a digital elevation model (DEM) with a pixel of 2 m derived from LiDAR-PNOA-cob2 2015 CC-BY 4.0 (scne.es) with a density of 0.5 point per m⁻².

Four of the largest saladas were selected as representative of the diverse wetland types recognized in the endorheic system (Castañeda et al., 2013): the playa lakes Rollico (40.8 ha) and Guallar (14.9 ha) with bare bottom, and the wet basins Gramenosa (29.9 ha) and Agustín (62.5 ha) with perennial halophytes and arable land. All them are affected by advancing plowing and cropping as the soil wetness allows the machinery to enter. Agustin occupies the bottom of a flat-bottomed valley with a relatively smooth topography and a higher agriculture imprint.

Two different methods were used for the soil salinity survey. A first campaign consisted in soil sampling along transects extending through the wetland-crop interface. The transects include the bottom (bare or vegetated), fringe of natural vegetation, and agricultural plots. According to the basin geomorphology, these transects were of NW-SE direction in Guallar and Rollico, and of NE-SW direction in Gramenosa (Figures 2a–c). The second soil sampling campaign was conducted in Agustin, in agricultural plots declared as “arable land” in the Land Parcel Identification System (LPIS), the information system recording all agricultural parcels in the EU Member States and key control mechanism under the Common Agricultural Policy (CAP) (European Court of Auditors, 2016). The accessible fringing area of the wetland was sampled with the EM sensor and three of these surveyed plots, named A1, A2, and A3, were selected for soil sampling (Figure 2d). Table 1 summarizes the soil sampling characteristics. A total of 34 sites were sampled along the 1870 m of the three transects and 28 sites were distributed in three plots of Agustín over a total of 34.4 ha.

Soil samples were extracted by using a 7 cm diameter Eijkelkamp hand auger at 20 or 25 cm intervals. At the different sites, the sampling depth and/or the depth intervals were not the same (Table 1) due to their different surveys or operators. The sampled depth of the soil ranged from 50 to 260 cm. A total of 183 soil samples were taken along the transects and 136 distributed in agricultural plots, totalizing 319 soil samples.

We measured the groundwater depth and used a bailer-type sampler for taking water samples. Electrical conductivity (EC) and pH were determined in the field with the ORION 150A + conductivity meter and CRISON PH25 pH meter, respectively.

Electromagnetic sensor (EMS) readings were taken with a hand-held Geonics EM38 sensor in points along wetland-crop transects in Rollico, Guallar, and Gramenosa. The readings were taken 9 days after a rain of 32.2 mm in the nearby Valfarta weather station (SIAR network ² ), in trafficable soil conditions. The rainfall was enough for the correct operation of EMS. The readings were taken on the same day in horizontal and vertical positions of the sensor, i.e., coils parallel and perpendicular to ground, respectively. The readings were corrected with the conversion table of US Salinity Laboratory Staff (1954) in order to reference them to 25 °C. Then these dimensionless numbers (Herrero et al., 2024), divided by 100, were termed EMH and EMV, respectively. In Agustin we used the electromagnetic sensor DualEm 1S operating in horizontal (EMh) and vertical (EMv) positions. After Urdanoz and Aragüés (2012), this sensor gives measurements basically similar to EM-38. The Dualem sensor was mounted on a plastic sled towed by a tractor at a distance of 3 m allowing the automatic acquisition of georeferenced readings. The sensor was connected to a Garmin Etrex GPS and an Allegro CX computer (Juniper Systems, USA), where the sensor readings were recorded and georeferenced in real time. A preliminary map of mobile EMS readings was created in the field using the ArcGIS^© kriging tool in order to decide the soil sampling points.

EMh y EMv were used for estimating soil salinity as ECe and EC1:5 at different sampled depth intervals using linear regression (López-Bruna and Herrero, 1996; Nogués et al., 2006) up to 100 cm. For the plots, the best regression equation of ECe was applied on the raster image produced with the DualEm 1S horizontal readings to generate the corresponding soil salinity map.

Soil samples were dried at 30 °C–40 °C to avoid the loss of water from the gypsum (CaSO₄•2H₂O) and grounded and sifted through a sieve of 2 mm mesh diameter. In all the soils we determined the electrical conductivity (EC) of 1:5 and 1:10 soil:water extracts with the conductivity meter Orion 013605MD. The pH was measured using a pH electrode (Orion 9157BNMD), and the major ions were determined by chromatography (Metrohm 861 Advanced compact IC, Metrohm AG, Herisau, Switzerland). The EC of the saturation extract (ECe) was determined in all the samples of the transects and in a half of the samples (69) of Agustin and it was estimated for the rest of Agustin samples using the regression of ECe on EC1:5.

The salinity profiles were analysed based on the EC1:5 extracts per sampling site using the original sampling depths. In order to compare the different sites, for each site we established a comparable depth up to 50 cm, as relevant for the cereal and vegetation roots activity, and to 100 cm to analyse possible influence of deeper layers or groundwater. The integrated salinity for these two depths was calculated by weighting the EC values by the thickness of each layer (Castañeda et al., 2012).

The detailed geomorphological maps of the saline wetlands show a diversity of landforms, including lacustrine deposits associated with the flooded or vegetated bottoms and aeolian deposits on the leeward side of some of the saladas (Figure 3). The four studied saladas are inset into broader depressions usually framed by conspicuous flat-bottomed valleys excavated on the flat Miocene structural surface. The valleys network evidences structural and hydrological connection between neighboring depressions which stand out for their elongated shape in a NW-SE direction. Guallar (Figure 3a) develops along a NW-SE elongated depression with a major axis of about 1 km long and delimited to the northeast by a gypsum escarpment of about 10.4 m height that includes decimetric limestone strata. The flat bare lakebed extends to the aeolian deposit mainly covered by halophytes.

Rollico (Figure 3b) is inset within a relatively deep depression developed in the confluence of a wide structural surface developed to the NW and a pediment slope draining from the SW. The bare lakebed is enclosed by two lacustrine terraces, both of them less than 6 m high above the lake bottom, and by two adjacent aeolian deposits in the leeward zone. Gramenosa (Figure 3c) is part of a wider depression which includes a continuous colluvial deposit surrounding the main wetland and a small adjacent lakebed at the NW zone surrounded by vegetation. Agustin (Figure 3d) is located at the lowest topographic position of a complex area excavated into the largest flat-bottomed valley in the region, about 22 km long and 1.7 km of maximum width. The border of the depression hosting the wetland is blurred within the valley which includes a noticeable colluvium with a landslide to the north, together with small pediments. Leeward aeolian deposits are also present in the bottom of the depression.

The scarce elevation differences accounted along the transects (Table 1) provide evidence of the flatness of the saladas bottoms, especially in the bare soil of the playa lakes. A relatively higher topographic gradient is present along the vegetated fringe, varying from +2.4 m in Guallar to +8.1 m in Rollico. The fringe of Agustin around the vegetated bottom is relatively flat, with less than 2.5 m of elevation difference between the sampled sites.

The regression of the EC1:5 versus EC 1:10 values obtained a R² = 0.98 in the transects and R² = 0.97 in Agustin plots. In both cases, the group of points more deviating are from line and near the origin of coordinates, corresponding to the low salinity values, less than 2.5 dS m⁻¹. A line can be fitted well (Figure 4) even without eliminating the points with lower EC, which probably show the presence of gypsum or maybe of more soluble salts in an insufficient proportion to reach the concentration of the 1:5 extract in the 1:10 extract.

In Agustin plots, the regression of ECe on CE1:5 for the 69 samples with both determinations gives a determination coefficient of 88% (Figure 5), and a standard error of 6.8 dS m⁻¹ and allows to estimate the ECe for all the Agustin samples with Equation 1: E C e = 8.0877   x   E C 1 : 5 − 11.948 (1)

Table 2 summarizes the soil salinity values for each salada as a whole and for each soil cover. The salinity shows a similar pattern in the three transects, in close relationship to the geomorphic position and the soil cover. The highest values occur in the bare soil bottom of the saladas (mean EC1:5 of about 15 dS m⁻¹), the intermediate values are found in the fringe of vegetation (from 6.5 to 10.2 dS m⁻¹), and the lowest values correspond to the soils of the crops (from 3 to 4.6 dS m⁻¹). Apart from the averaged values, there are also differences in the range of the EC values along the transect in each salada. The transects along Guallar and Rollico playa-lakes collect the highest salinity range (up to 27.7 dS m⁻¹) and a similar mean salinity at their bare bottom. The lowest mean salinity occurs in the vegetated saladas, Gramenosa and Agustin, with EC1:5 of 5.7 and 4.6 dS m⁻¹, respectively. Here the salinity range is about a third part of the range in the playa-lakes, about 9 dS m⁻¹.

The natural vegetation of Rollico fringe (Table 2) withstands higher salinity than that of Guallar or Gramenosa, not only by the mean but also the minimum EC values registered. The crops sampled along the transects exhibit also notable soil salinity, especially high in Agustin (mean 4.6 dS m⁻¹), followed by Guallar and Gramenosa, with 3.6 dS m⁻¹ and 3.1 dS m⁻¹, respectively.

Based on the ECe salinity phases of Nogués et al. (2006), a 73% of the soil samples are very strongly saline (ECe >16 dS m⁻¹). They are spread in the saladas regardless of the soil cover, from the hypersaline bare soil of the playa-lakes to the crops. Similarly, most soil samples of the natural vegetation are very strongly saline. All the soil salinity phases, except the non saline, can be found in the cereal area: slightly (5%), moderately (14%), strongly (23%) and very strongly (58%) saline. Summarizing, half of the very strongly saline soils samples have been taken in agricultural areas with no soils free of salinity.

Figure 6 illustrates the vertical distribution of soil salinity (EC1:5) within each salada, along the transects and in the dryfarmed area fringing Agustín. Salinity profiles were analyzed to depths ranging from 60 to 260 cm. The soil depth reached was very variable and, in general, lower at sites located on the crop area of the transects. At Guallar (Figure 6a) the sampled depth uniformly decreased towards the crops whereas at Rollico (Figure 6b) a similar depth was sampled in all the sites along the transect. There is a noticeable variability of the vertical distribution of soil salinity along the transects though in general, the salinity profiles can be grouped by soil cover, especially the upper soil horizons. Deeper horizons show a heterogeneous salinity pattern since the changes may occur at different depths. In general, the range of the soil salinity values along each transect is maximum at the surface horizons and clearly decreases with soil depth, especially in Guallar and Rollico.

The salinity profiles of Agustín reached up to 150 cm having a relatively uniform sampled depth (Figures 6d–f). The vertical distribution of salinity is heterogeneous even in the same plot and presents no relationship to the elevation. Unlike the transects, the range of salinity values are similar at the surface and subsurface horizons. The highest depth (150 cm) and the highest salinity range, EC1:5 from 2.5 to 11.5 dS m⁻¹, was reached at A1 plot (Figure 6d). The salinity can be maximum at the surface horizon, constant along the profile or can increase with depth. Frequently, a higher salinity ranging from 3.4 to 9.4 dS m⁻¹ is observed in the first 50 cm. Only one point showed a considerable salinity increase deeper than 100 cm.

The saline profiles of A2 were quite uniform and represent a lower salinity range, from 2.3 to 6.6 dS m⁻¹, with only two sites with EC increasing with depth. All the A3 plot profiles were slightly shorter (up to 120 cm) and, similar to A1, showed a diversity of vertical salt distribution and the surface horizon (25 cm) can be very strongly saline, with EC1:5 up to 7.3 dS m⁻¹.

Groundwater was reached in some of the sampled sites along the transects (Figure 7), at depths ranging from 31 to 195 cm (before equilibrium) plus one site flooded in Rollico bottom. Once stabilize, the water level was at ∼140 cm in sites #5 to #9 in the vegetated bottom of Gamenosa, with a mean EC of 82.6 dS m⁻¹; up to 25 cm at the bare bottom of Guallar (sites #1 to #3) and deeper (∼190 cm) in the vegetated fringe (site #4), with 133 dS m⁻¹ mean EC; and at ∼70 cm in the bare bottom and vegetated fringe of Rollico (sites #1 to #6), with a mean EC of 144.3 dS m⁻¹. In general pH ranged from 6.8 to 7.3 and the anions composition was characterized by the high amount of chloride (54 g L⁻¹) followed by sulfate (27 g L⁻¹) with a ratio Cl/SO₄ of 2.7. The main cations were sodium (21 g L⁻¹) followed by magnesium (8 g L⁻¹) with a mean ratio Mg/Ca of 20.1.

Figure 8 summarizes the soil salinity of the saturation extract (ECe) weighted up to 50 cm (ECe50) and 100 cm (ECe100). The overall distribution of EC values at the sampled sites was consistent with the expected salinity. A limited number of sites can be classified as slightly or moderately saline. The salinity range along the transects showed some differences between both soil depths. In general, ECe100 was higher than ECe50, and the EC range of the vegetated saladas was not as wide as that of the playa-lakes. The ECe50 of the playa-lakes bottom can be twelve times that of the crops in Guallar and more than twenty times in Rollico. In these playa-lakes the ECe50 gradient was greater than that of ECe100 due to the high ECe50 values at their hypersaline bottom. The EC variability within Agustin plots can not be associated with systematic relevant topographic differences or the phreatic level.

The graphs of Figure 8 show some local spots of unexpected high salinity, identified by the ECe of both depths. Most of these deviated points are located outside the salada bottom, in the leeward low-lying area of the aeolian deposits (orange in Figure 3) identified in Guallar (Figure 8a, site #7), Rollico (Figure 8b, sites #9, #10), and Agustín A2 (Figure 8e, sites #4 to #6). A similar saline spot occurs on the edge at Gramenosa (site #4) (Figure 2). The water table was never reached at these sites during sampling or during the visits in the following 3 months.

The calibration of EMS to determine edaphic salinity in the studied transects, using the EC obtained from 1:5 and 1:10 soil extracts for integrated soil depths resulted in the equations summarized in Paracuellos (2006). The correlation analysis showed differences between the saladas (Figure 9). Only Rollico (Figures 9b,e) did not show differences in the coefficients of determination (R²) between horizontal and vertical readings. At Gramenosa (Figures 9c,f), R² values were very low, the maximum being 4% in extract 1:5 at a depth of 0–20 cm. In both extracts, the regressions with EMh had higher determination coefficients than those with EMv for all depths. In 1:5 extract the best fit occurred up to the depth of 20 cm and became worse as depth increased, whereas in 1:10 extract the best fit occurred for the depth of 60 cm and thereafter R² decreased with soil depth. At Guallar (Figures 9a,d), R² values were acceptable for soil depths down to 80 cm (74%–87%), although there was a sharp decrease from this depth (45%–63%). The R² values were higher for settings with EMv (∼81%) than with EMh (∼69%). In Guallar the prediction EC1:5 values were established with the equation Equation 2, E C 1 : 5 = 2.93 + 1.07 × E M v (2)

At Rollico (Figures 9b,e), for both extracts R² increased with depth even though the number of points that influence is smaller and the standard error is also smaller (Paracuellos, 2006). The maximum R² occurred at 100 cm, 94% and 96% at 1:5 and 1:10, respectively. R² was higher at all depths in the 1:10 extract. No significant differences in R² were observed in EMh and EMv.

For Guallar and Rollico the calibration equations corresponding to the 60 cm depth were appropriate because the salinity of this thickness was the most influential on cereal yield. They showed acceptable R², especially the 1:5 extracts (Figures 9a,b), whose values were 80% or higher. In Rollico, the calibration equation is Equation 3, E C 1 : 5 = 2.53 + 0.63 × E M h (3)

The map of EMS readings (Figure 10) displayed a noticeable spatial variability within the Agustin plots. Higher R² values were obtained for the EMh than for EMv. The R² values were very low at the first layer (0–25), increasing from 0 to 50 cm. The correlation with EC1:5 was slightly better than with EC1:10. We selected the equation Equation 4, E C e = 3.095 + 6.561 × E M h (4) with R² = 81%, corresponding to the weighted ECe up to 50 cm depth, the zone where most roots responsible for nutrient absorption under irrigation usually grow (Isla, 2001).

The resulting salinity map provided soil salinity estimation at unsampled sites (Figure 11a) and allowed for a comprehensive description of the salinity distribution and comparisons between the surveyed areas. Almost the complete fringe of Agustín depicted ECe >4 dS m⁻¹. The salinity decreased towards the crops in areas with higher elevation, especially in A2, because at higher topographic positions the depth of the water table increased. This area exhibited the greatest spatial variability of salinity from the studied plots. A concentric distribution was observed, with ECe >10 dS m⁻¹ at the center and the rest of the area with values between 2 and 10 dS m⁻¹. The higher values occurred at the leeward depression of the aeolian deposit (Figure 3d) currently cultivated though with limited production (Figure 11b).

A diversity of landforms was identified in the saladas area, covering about 230 km², including halophytic vegetation in the saladas bottoms, aeolian deposits associated with bare playas, sharp escarpments, nested or coalescent depressions and networks of flat-bottomed valleys, among others. The detailed geomorphological map, based on the aerial photographs from 1957, captured the geoforms under landscape conditions previous to the homogenization produced by the transformations carried out to install the new irrigation (Domínguez-Beisiegel et al., 2013a).

A smoothed and homogeneous landscape was found during the soil surveys especially due to the land consolidation and excavation works for irrigation and drainage pipes. Following the disappearance of many of the saladas, the area now shows larger and more continuous agricultural fields free of vegetated fringes and with a great number of new intensive pig farms (Tierra et al., 2025). Under this scenario, the original mapped landforms (Figure 3) are key for interpreting the distribution of soil salinity and the frequent occurrence of superficial color patches with distinct soil properties (Castañeda and Moret-Fernández, 2013).

Similar transformations occurred in other semiarid environments in the CEB, such as the saladas of Alcañiz-Calanda, among others, and in other regions of Spain and saline environments worldwide usually requiring in situ knowledge to predict salinization effects and prevent degradation (Tedeschi, 2020).

The EC of the saturation extracts is a worldwide accepted method established by the US Salinity Laboratory Staff (1954), originally aimed at agriculture for expressing soil salinity. We applied it in wetland ecosystems for meaningful comparisons in the continuum of habitats at the interface with agriculture. The EC determined in the two extracts, 1:5 and 1:10, are robust as showed by the relationships between them. The relationship between the ECe and EC1:5 is similar to that obtained by Herrero et al. (2015) computing 359 samples from 10 saladas not coincident with the ones studied in the present article. The regression result of ECe on EC1:5 in Agustín plots (R² = 88%) is comparable to that obtained for soils from the vegetated bottom of Agustín by Herrero (2008).

The response of EMS in Agustin presents the best relationship with a soil depth of 0–100 cm in A3 and with EMh readings (Figure 9f) although the relative contribution of the deep soil layers to the SEM readings is greater for the vertical position of the receiving coil (Abdu et al., 2007). These authors established a sensitive depth of exploration twice for vertical than for horizontal modes in a similar sensor. Yao and Yang (2010) proposed independent equations for each depth but Huang et al. (2015) considered that it is not necessary.

In Gramenosa, no acceptable fit was obtained, with low R² values. Despite this, the graphs obtained maintain a relationship in their form, and the regression equations show a certain consistency, indicating that the results could be more favorable if the number of observations increased. In Guallar, the fit was generally good, with R² values sharply dropping from 80 cm onward, possibly due to the smaller number of points or the influence of deeper layers. This is confirmed by the higher R² values in EMv, where the relative contribution of deeper layers is greater. Rollico provided a good fit for all depths and with great uniformity, even though the number of computed points decreased with depth.

The transect method is an appropriate approach to explore the soil salinity in the flat landscape of Monegros with scattered depressions. The inclusion of the gradients of salt-tolerant vegetation in the transects made it possible to capture a representative range of salinity. The resulting salinity range detected along the studied transects (from 9 to 26 dS m⁻¹) (Figure 6) can be used as a model to be extrapolated to all the area hosting about 150 depressions and halophytes (Conesa et al., 2011; Castañeda et al., 2013) with the help of digital mapping techniques together with soil salinity proxies (Zhang et al., 2011).

The geomorphological features, including topography, govern the movement and distribution of soil salinity as they condition how and where water flows accumulates. In this sense, the spatial distribution of salinity together with appropriate geomorphological maps (Figure 3) are very useful to predict the effects of irrigation and salt mobilization at the landscape scale (Duncan et al., 2008). On the other hand, the salinity maps based on calibrated mobile EMS readings are advantageous for continuous information of soil salinity, though the vehicle approach is restricted to accessible areas, e.g., agricultural plots or other semi-urban environments (Williams et al., 2024), not frequently found in natural areas.

The integrated salinity (Figure 8) allows comparing sites, soil covers, and saladas. The irregular soil salinity in relatively small agricultural plots (Figures 8d–f) with no relevant elevation differences could be related to the proximity of deeper saline layers even if no groundwater was reached on the sampling dates. García-Vera (1996) identified in Agustin the confluence of two aquifers in the area having dissimilar characteristics. Salvany et al. (1996) documented the frequent change of geological materials and groundwater characteristics at distances of few meters in the area, and the occurrence of local increases in salinity at the saladas fringe due to variable-density effects of groundwater. Besides, a great number of tectonic joints and the occurrence of preferential flows probably contributes to the heterogeneity of the subsurface materials yielding differences in soil salinity at distances of just a few meters.

Regarding the salinity profiles of Figure 6, one of the differences between the playa lakes and the agricultural plots is the striking high salinity range of the upper horizons in natural areas. The high ET₀ favors the capillary rise and evaporation through the bare soil (Castañeda and García-Vera, 2008) and subsequent salt accumulation in the upper layers, generating hypersaline areas that intensify the salinity gradient. On the contrary, the agricultural practices in the plots of Agustin very probably have a homogenizing effect on the upper layers of the soil profile (Ohigashi et al., 2025). To comply with CAP measures, continuous tillage is required under rainfed, even under long fallow practices (16–18 months).

Low values of Gramenosa may depend on the points chosen. Irregular salinity values in vegetated areas may be related to the higher density of the vegetation cover and the abundance of nitrophilous plants interspersed with the salinity-driven plant succession. Under these conditions, the soil cover changes cannot be precisely distinguished when selecting sampling sites. Sampling more closely aligned with lateral variation in soil salinity would require a prior study of plant distribution at a very detailed scale.

The salinity intervals for soil salinity classes (Nogués et al., 2006) were designed for agricultural soils and leave the samples with ECe >16 dS m⁻¹ in a single group, that is, 73% of the 319 samples studied. Regarding these intervals, new considerations would be needed to integrate the agricultural and ecological perspectives into the study of soil salinity in hypersaline areas.

The surface extent affected by salinity of Figure 11a can be quantified for individual CAP parcels of the interface between crops and the protected habitats. For agronomic application, a production criterion could be applied to establish a reference limit of production for each parcel (Serrate, 2009). As an example, Royo and Abió (2003) established as criterion a 50% decrease in durum wheat and barley yield due to the presence of salts in the soil. The authors obtained a mean salinity value (EC50) of 11.3 dS m⁻¹ as a limiting factor. EC50 is the electrical conductivity of saturation extract, which reduces yield to half of what would be obtained under non-saline soil conditions. This value would be more restrictive than that of barley, which can tolerate much higher salinity levels, with EC50 = 18 dS m⁻¹ (Ayers and Westcott, 1985). In this context, a proposal for a new agri-environment-climate measure could include qualifying plots with ECe >10 dS m⁻¹ in more than 50% of the surface area. These plots could be considered not suitable for cultivation such as that of Figure 10b.

Some of the CAP greening obligations are to be monitored using the Land Parcels System (LPIS). Since part of the CAP subsidies to farmers are to be paid for agricultural practices beneficial for the climate and environment (e.g., benefit biodiversity, habitat restoration and maintenance, or retain landscape features), the plots of limited agricultural potential may form part of a new category of eligible areas within the LPIS as non-productive areas prone to natural variable salinity as proposed by Serrate (2009). Moreover, as recommended by European Court of Auditors (European Court of Auditors, 2016) to accurately identify this type of land as eligible, targeted on-the-spot data may be required, and this requirement is well suited for soil salinity due to their intrinsic temporal variability.