Impact of Irrigation Management Decisions on the Water Footprint of Processing Tomatoes in Southern Spain

Abstract

The water footprint is an increasingly demanded environmental sustainability indicator for certifications and labels in agricultural production. Processing tomatoes are highly water-intensive, and existing studies on water footprint have uncertainties and do not consider the impact of different irrigation configurations (e.g., surface drip irrigation (SDI) and subsurface drip irrigation (SSDI)) and irrigation strategies. This study presents a two-year experimental investigation to determine the water footprint of processing tomatoes grown in southern Spain (Andalusia) and the impact of SSDI and deficit irrigation. Five irrigation treatments were established: SDI₁ (surface drip irrigation without water limitation), SDI₂ (surface drip irrigation without water limitation adjusted by soil moisture readings), SSDI₁ (subsurface drip irrigation without water limitation and a dripline depth of 15 cm), SSDI₂ (similar to SSDI₁ but with mild/moderate water deficit during the fruit ripening stage), and SSDI₃ (subsurface drip irrigation without water limitation and a dripline depth of 35 cm (first year) and 25 cm (second year)). Measurements included crop vegetative growth, leaf water potential, leaf gas exchange, nitrate concentration in soil solution, and crop yield and quality. The soil water balance components (actual evaporation, actual transpiration, deep drainage), necessary for determining the total crop water footprint, were simulated on a daily scale using Hydrus 2D software. Results indicated that SSDI makes more efficient use of irrigation water than SDI. The water footprint of SSDI₁ was 20–35% lower than that of SDI₁. SSDI₂ showed similar water footprint values to SDI₁ under highly demanding environmental conditions and significantly lower values (≈40%) in a year with lower evaporative demand. The dripline depth in SSDI was critical to the water footprint. With a 35 cm installation depth, SSDI₃ had a significantly higher water footprint than the other treatments, while the values were similar to SSDI₁ when the depth was reduced to 25 cm.

1. Introduction

The concept of ‘water footprint (WF)’ has been defined to raise social awareness of water usage and to analyze the relationship between human consumption habits and their impact on natural resources. It serves as an indicator of water resource utilization in human activities by assessing both the quantity of water used and its impact on water quality [1,2,3,4]. The WF of a product reflects the volume of freshwater used to produce a specific product, measured throughout the entire supply chain [5]. Since most of the water usage for many food products arises during the agricultural production stage [6], recent endeavors have focused on quantifying the WF of agricultural practices to use it as a decision-making tool concerning the utilization and allocation of water resources in stressed areas [7].

The total WF of agricultural crops can be calculated as the sum of three components, i.e., blue WF, green WF, and grey WF [5]. Blue WF refers to the consumptive use (mainly evaporative) of water from surface water and groundwater sources. In irrigated agriculture, this component would reflect the amount of water applied through irrigation that has been used in the evapotranspiration process. Blue water volumes do not necessarily coincide with the volumes of applied irrigation, as the term ‘water use’ refers only to water used in the evapotranspiration process and that does not return to its source [5]. On the other hand, green WF refers to the amount of rainwater stored in the soil that is used in the evapotranspiration process, while grey WF refers to the volume of freshwater needed to dilute pollutants discharged into water bodies because of irrigation.

Processing tomato (Solanum lycopersicum L.) is a high-value crop in Spain, which is the fourth-largest producer globally, following the United States, China, and Italy, with an annual average production of nearly 3 million tons, representing around 7–8% of global production [8]. Andalusia (Spain) contributes to approximately 30% of the national processing tomato production, covering a cultivated area of approximately 6000–7000 hectares [9]. Processing tomato is a highly water-demanding crop, with water requirements that typically exceed 600 mm per growing cycle in Spain’s main producing areas [10]. Due to its high water consumption, in drought-prone areas like Andalusia, processing tomatoes are among the first crops to suffer the impact of water restrictions. Thus, there is significant interest within the sector to establish production and crop management systems that optimize water usage and reduce the crop’s WF. However, recent studies on the water footprint of processing tomatoes in Andalusia are lacking, underscoring the importance of measuring and evaluating the WF to ensure sustainable water management in this agricultural region.

It is known that the WF of crops is conditioned by factors such as the region’s climate, production system, and soil type. While in Morocco the WF of one kilogram of fresh tomato production in unheated greenhouses does not exceed 30 L [11], in Italy, the production of one kilogram of processing tomatoes uses, on average, 114 L of water [3]. A former study conducted in Spain determined that the WF of producing one kilogram of tomatoes varies between 216 and 306 L, with an average value of 236 L kg⁻¹ [12]. However, these figures are much higher than those observed in Greece, where processing tomato production has a WF ranging between 37 and 131 L kg⁻¹, depending on the climatic conditions and soil characteristics of each area [13].

Despite having estimations of the WF of processing tomatoes in various Mediterranean regions, the high variability in the figures generates some uncertainty, probably due to the lack of reliability or absence of robust data on which to calculate the WF of agricultural production [14]. A clear example is the fact that less than half of the studies on WF include the determination of grey WF, and those that do rely on simplifications or estimates (e.g., assumption of 10% of nitrogen fertilizer losses through drainage or runoff) rather than measured or estimated data in experimental plots. Another weakness is that most WF studies estimate crop evapotranspiration (ET) under standard conditions (i.e., non-limiting soil water conditions and surface drip irrigation), without analyzing the impact that other irrigation scenarios (e.g., deficit irrigation and/or subsurface drip irrigation) may have on ET and hence the crop’s WF.

Both subsurface drip irrigation and deficit irrigation strategies have shown great potential as a tool for reducing irrigation inputs with little or no impact on the production of certain crops [15,16]. Deficit irrigation involves supplying water below the crop’s requirements, either throughout the entire growing season (sustained deficit irrigation) or during specific, less critical phenological stages (regulated deficit irrigation). In the context of processing tomatoes, deficit irrigation has not always shown an increase in water use efficiency compared to full irrigation [17,18], although when mild to moderate water deficits were applied only during specific phenological periods, such as from the early stages of fruit growth onward, the performance of deficit irrigation surpassed that of full irrigation [10,19].

Subsurface drip irrigation (SSDI) has proven to be an irrigation method with high potential to improve water use efficiency on several crops. The reduction of the evaporative component of the soil together with a higher wetted soil volume and better use of water and nutritional resources seem to be the reason for this behavior. In the review conducted by Lamm [20], it is observed that SSDI increased the production of processing tomatoes with respect to surface drip irrigation (SDI) by 7%, on average. However, they also highlighted significant yield reductions in light-textured soils and the need to supplement SSDI irrigation with a sprinkler irrigation system to assist in the early stages of crop establishment due to the usual depth at which pipes are buried in these systems. When deficit irrigation and SSDI are combined, water use efficiency can be even greater compared to full irrigation either with or without SSDI [21].

Based on this evidence, it is necessary to design experimental trials that allow elucidating the impact that these irrigation management strategies can have on the WF of processing tomatoes. This study has proposed an experimental investigation to evaluate the agronomic performance of processing tomatoes under different drip irrigation configurations (surface vs. subsurface drip irrigation) and management strategies (full irrigation vs. deficit irrigation). Crop development data as well as the actual irrigation applications have been used to estimate the components of the soil water balance daily, which are necessary for estimating the WF of tomato cultivation under the indicated management scenarios and under the edaphoclimatic conditions of the Middle Guadalquivir River Valley.

2. Materials and Methods

2.1. Description of the Experimental Field

The research was conducted on a commercial processing tomato farm located in Southern Spain (Carmona, Andalusia; 37.505978° N, −5.845658° W) over two growing seasons (2022 and 2023). Processing tomatoes (cv. Heinz 1015) were grown on a uniform sandy clay loam soil (Sand = 60%, Silt = 17%, Clay = 23%), with initial average values of pH 7.8, electrical conductivity 158 μS/cm, and organic matter content 0.6%. The crop was planted on ridges 1.10 m wide and 0.2–0.3 m high, with a spacing of 1.5 m between centers. Transplanting occurred on 18 April 2022 and 3 April 2023. A single planting line per ridge was set up, with a planting density of approximately 4 plants per meter of ridge length. The climate in the area is Mediterranean, with rainfall primarily concentrated between October and May, and an average annual precipitation and reference evapotranspiration (ET₀) of 496 ± 182 mm and 1395 ± 73 mm, respectively, based on 2001–2023 averages recorded at a nearby agroclimatic station managed by the Andalusian Government.

Table 1. Climatic variables registered at the experimental site during the 2022 and 2023 growth cycles. DAT: days after transplanting; T_m: daily mean temperature (°C); RH_m: daily mean relative humidity (%); u: wind speed at a height of 2 m (m s⁻¹); R_s: solar radiation (MJ m⁻² d⁻¹); P: rainfall (mm); ET₀: reference evapotranspiration calculated with the Penman–Monteith equation [22] (mm).

Growth Cycle 2022
DAT	Tm	RHm	u	Rs	P	ET0
0–10	15.8	59.6	2.3	20.5	35.8	47.4
11–20	19.6	55.4	1.5	26.0	21.6	52.5
21–30	23.5	45.7	1.6	27.7	0.0	64.1
31–40	23.9	42.4	1.9	27.7	0.0	69.7
41–50	22.4	51.4	2.1	27.4	0.0	63.0
51–60	29.0	38.6	1.5	27.9	0.0	74.3
61–70	22.7	44.9	2.4	28.8	0.0	70.5
71–80	25.2	43.9	1.6	29.4	0.0	69.7
81–90	30.8	32.0	1.6	28.8	0.0	82.2
91–102	29.5	37.6	1.3	28.3	0.0	84.4
Growth Cycle 2023
0–10	18.6	39.2	1.1	21.4	0.0	45.0
11–20	19.7	42.9	1.6	26.1	1.6	55.5
21–30	23.1	46.5	0.8	25.4	0.0	50.3
31–40	22.5	45.5	1.0	27.4	0.0	55.4
41–50	18.9	54.0	1.2	22.9	19.6	46.8
51–60	19.7	67.4	0.9	21.2	31.6	41.8
61–70	22.1	64.4	1.1	22.1	47.2	46.6
71–80	25.2	51.9	1.3	28.1	0.0	61.7
81–90	28.7	43.0	0.9	29.0	0.0	65.1
91–100	28.0	45.6	1.3	29.0	0.0	69.1
101–105	29.3	38.0	1.2	29.8	0.0	42.2

shows the climatic data recorded during the study period at an agroclimatic station belonging to the network of agroclimatic stations of the Junta de Andalucía, located near the experimental plot. In 2022, the crop cycle lasted 102 days, with average records of temperature, relative humidity, wind speed, and solar radiation of 24.2 °C, 45.1%, 1.8 m s⁻¹, and 27.3 MJ m⁻² d⁻¹, respectively. The crop cycle was characterized by three rainfall episodes in the first 20 days, providing a total of 57.4 mm of water to the crop, and a thermal regime marked by two heatwaves in June (days after transplanting, DAT, 51–60) and July (DAT 81–102), with average daily maximum temperatures of 38.5 °C and 39.9 °C, respectively. The accumulated reference evapotranspiration (ET₀) throughout the 2022 crop cycle was 678 mm. In 2023, the crop cycle lasted 105 days, with average records of temperature, relative humidity, wind speed, and solar radiation of 23.3 °C, 44.9%, 1.1 m/s, and 25.7 MJ m⁻² d⁻¹, respectively. The crop cycle was characterized by several rainfall episodes mid-cycle, providing a total of 100 mm of water to the crop. The accumulated reference evapotranspiration (ET₀) during the 2023 crop cycle was 580 mm, a 15% decrease compared to the values observed in 2022.

The crop was irrigated with 16 mm diameter polyethylene driplines delivering a flow rate of 1.6 L h⁻¹ per emitter, spaced 0.5 m. The drippers, model AZUD GENIUN (SISTEMA AZUD S.A., Alcantarilla, Spain), were pressure-compensating and anti-suction. A single dripline was established for each row of plants. The irrigation head consisted of the pumping system, a self-cleaning disc filtration system model AZUD HELIX AUTOMATIC, and a fertigation equipment model AZUD QGROW (SISTEMA AZUD S.A., Alcantarilla, Spain).

2.2. Irrigation Treatments

Five irrigation treatments with four replications per treatment, each covering approximately 30 × 40 m², were established. Among the five irrigation treatments, two were surface drip irrigation (SDI) treatments (SDI₁ and SDI₂) and three were subsurface drip irrigation (SSDI) treatments (SSDI₁, SSDI₂, and SSDI₃). Treatments SDI₁ and SDI₂ were irrigated to meet the crop’s water needs (ETc), with the exception that the irrigation doses in SDI₂ were adjusted based on soil moisture readings recorded with a 90-cm-deep profile probe model Drill and Drop (Sentek Sensor Technologies, Stepney, Australia). The probe has sensors every 10 cm, from a depth of 5 cm down to 85 cm. Of the three SSDI treatments, two had the driplines installed at a shallower depth (SSDI₁, SSDI₂) and one at a greater depth (SSDI₃). The dripline depth for SSDI₁ and SSDI₂ was 15 cm, whereas depths of 35 cm in 2022 and 25 cm in 2023 were used in SSDI₃. This variation between years was made to reduce the percolation losses observed in SSDI₃ in 2022.

Regarding the water doses applied, SSDI₁ received 90% (85% in 2023) of the water applied to SDI₁ throughout the entire crop cycle. SSDI₂ was scheduled like SSDI₁ except during the fruit maturation phase (the last month of the crop cycle) where it received 80% (70% in 2023) of the water applied to SDI₁. SSDI₃ received the same amount of water as SDI₁ in 2022 and 85% of SDI₁ in 2023.

Table 2. Description of the irrigation treatments. SDI: surface drip irrigation; SSDI: subsurface drip irrigation.

Treatment	Type	Dripline Depth (cm)		Irrigation (% ETc)
		2022	2023	2022	2023
SDI1	SDI	0	0	100	100
SDI2	SDI	0	0	100 1	100 1
SSDI1	SSDI	15	15	90	85
SSDI2	SSDI	15	15	90|80 2	85|70 2
SSDI3	SSDI	35	25	100	85

¹ Soil-moisture-based adjustments of irrigation depths were performed. ² Percentages of ETc replacement during the vegetative and maturation growth stages, respectively.

summarizes the characteristics of the five evaluated irrigation treatments. Each irrigation treatment was equipped with a flow meter to record the daily volume of water applied. The flow meters were connected to a Redin model datalogger (Nutricontrol, Cartagena, Spain), which transmitted the data remotely on an hourly basis. Except during rainy periods, irrigation frequency was daily or even higher (2 irrigations per day) on days with high evaporative demand.

Crop water requirements were calculated using the FAO-56 approach [22] and the crop coefficients (Kc) derived for local conditions by Campillo-Torres et al. [23], which must be determined from the following expression:

$$K_c=-0.0001·{SS\left(\mathit{\%}\right)}^2+0.0202·SS\left(\mathit{\%}\right)+0.3033$$

(1)

where SS represents the area shaded by the crop. To determine SS, during the first two months of the crop cycle, coinciding with the phases of transplanting, flowering, and fruit set, and the early stages of fruit growth, weekly measurements of the fraction of soil covered by vegetation were conducted in all irrigation treatments. On sampling days, one RGB image per replicate was taken at selected locations using a mobile phone camera oriented toward the ground. An aluminum frame, specially constructed with internal dimensions of 1.50 m × 0.67 m, equivalent to an area of one square meter, was used to frame the RGB images (Figure 1a). The images were then cut to match the internal dimensions of the aluminum frame (Figure 1b) and subsequently processed using the Canopeo^® application (Oklahoma State University, Department of Plant and Soil Sciences, https://canopeoapp.com/) to determine the percentage of soil covered by the crop (SS, %) (Figure 1c). The Canopeo^® application has a resolution of 0.01% and an accuracy in the identification of pixels corresponding to green vegetation that exceeds 90% [24].

2.3. Nitrate Concentration in the Soil Solution

Although fertilization was the same for all treatments, a Rhizon soil moisture sampling probe (Royal Eijkelkamp, Giesbeek, The Netherlands) was installed per treatment at a root depth (30 cm) to estimate the average nitrate concentration in the root zone (Figure 2). This measurement is necessary for calculating the gray water footprint, as explained in the sections following. Soil solution was extracted to determine nitrate concentration using an RQFlex reflectometer (Merck KGaA, Darmstadt, Germany) on two dates in 2022 and three dates in 2023 during the crop cycle.

2.4. Crop Water Status and Leaf Gas Exchange Measurements

Crop water status was determined by measuring midday leaf water potential (LWP). These measurements were conducted using a Scholander-type pressure chamber (Figure 2), model SF-Pres-35 (Solfranc Tecnologías S.L., Tarragona, Spain), on one leaf per replication. The Water Stress Integral (WSI) was calculated using the following equation [25]:

$$WSI=\left|{\sum }_{i=0}^{i=t}\left({\overline{\gamma }}_{i,i+1}-c\right)n\right|$$

(2)

where t is the number of leaf water potential measurements taken, ${\overline{\mathit{\gamma }}}_{\mathit{i},\mathit{i}+\mathbf{1}}$ is the average value of the water potential measurements taken in the interval [i, i + 1], c is the maximum water potential observed during the measurement period (−0.55 MPa in 2022 and −0.47 MPa in 2023), and n is the number of days in the interval [i, i + 1].

Leaf photosynthesis rate (A) and stomatal conductance (g_s) were also determined in one leaf per replication using a portable gas exchange meter model CIRAS-3 (PP Systems International Inc., Amesbury, MA, USA) (Figure 2). For the gas exchange measurements, a photosynthetically active radiation flux density (PPFD) of 1500 μmol m⁻² s⁻¹ provided by an internal LED light source was established. The CO₂ concentration in the chamber was set at 400 μmol mol⁻¹. The gas exchange measurements could only be carried out during the 2022 season due to a failure of the equipment.

2.5. Yield

The harvest took place on 29 July 2022 (102 DAT), and 17 July 2023 (105 DAT). In each replication, a two-meter-long bed section was labeled and manually harvested (Figure 2). All tomatoes of commercial quality for processing were collected and weighed using a portable scale. After weighing, a random sample of approximately 2 kg of fruit per replication was taken, placed in zip-lock plastic bags, and kept in portable coolers until arrival at the laboratory, where they were stored in a refrigerator at 4 °C. In the laboratory, the tomato samples were blended with a mixer, and the total soluble solids content (°Brix) of each sample was determined using a portable meter (PAL-BX/Acid, Atago, Tokyo, Japan).

2.6. Soil Water Balance Components

The software Hydrus 2D/3D (v. 1.11) [26] was used to simulate the daily soil water balance components (actual evaporation, actual transpiration, deep drainage). The hydraulic properties of the soil were characterized using the van Genuchten–Mualem constitutive relations [27], for which the Rosetta model included in Hydrus was employed.

Table 3. Soil hydraulic parameters used in the modeling study. θ_r: residual water content; θ_s: saturated water content; α, n, and l: van Genuchten shape parameters; K_s: saturated hydraulic conductivity.

Soil Texture	θr (cm3 cm−3)	θs (cm3 cm−3)	α (cm−1)	n	Ks (cm d−1)	l
Sandy clay loam	0.0637	0.3865	0.0264	1.3321	15.55	0.5

presents the hydraulic parameters used in the modeling study assuming soil homogeneity within the flow domain (1 m depth).

Crop root distribution was defined using the model by Vrugt et al. [28]:

$$b\left(x,z\right)=\left(1-{\displaystyle \frac{z}{Z_m}}\right)\left(1-{\displaystyle \frac{x}{X_m}}\right)e^{-\left({\displaystyle \frac{p_z}{Z_m}}\left|z^{\ast }-z\right|+{\displaystyle \frac{p_r}{X_m}}\left|x^{\ast }-x\right|\right)}$$

(3)

where X_m and Z_m are the maximum rooting lengths in the x and z directions [L], respectively; x and z are the distances from the planting line in the x and z directions [L], respectively; px [-] and pz [-] are empirical parameters; x* [L] and z* [L] describe the location of maximum root water uptake in the vertical (z*) and horizontal (x*) directions and b(x,z) denotes the two-dimensional spatial distribution of potential root water uptake [-].

Table 4. Parameters of the root distribution model by Vrugt et al. used in the simulations.

Treatment	Zm (cm)	z* (cm)	pz	Xm (cm)	x* (cm)	px
SDI1, SDI2	45	15	0	45	15	0
SSDI1, SSDI2	45	15	3	45	15	0
SSDI3	45	25	3	45	15	0

and Figure 3 show, respectively, the root distribution parameters and root distribution patterns defined for the five irrigation treatments. The root distribution parameters used in the model are constant throughout the entire crop cycle. This approach is standard practice in soil water dynamic simulations for annual crops, particularly when using two-dimensional models [29].

The reduction in root water uptake under water stress was simulated using the equation by Feddes et al. [30], using the parameterization for tomato crops available in the database included in the Hydrus program. Due to the symmetry on both sides of the irrigation line, only half of the flow domain was simulated (Figure 3). The flow domain was defined as a rectangle measuring 75 cm by 100 cm, with a trapezoidal section (short base = 50 cm; long base = 60 cm; height = 20 cm) on the left side, representing the dimensions of the planting bed (Figure 3). This domain corresponds to a vertical plane perpendicular to the irrigation lines, bounded by the position of the planting line (coinciding with the irrigation line) and the midpoint between adjacent planting beds. The assumption of a two-dimensional (vertical) flow to model soil water distribution for drip irrigation systems is widely accepted [31,32,33] when emitter spacing is short and wetting bulb overlap is high, conditions that are met in the experimental setup of this study. The flow domain was discretized into a mesh of triangular elements, with refinement at the emitter location to ensure smaller elements in areas with rapid flow changes, allowing accurate representation in the numerical solution.

Regarding the boundary conditions defined at the outer edges of the flow domain, a no-flow condition was established on the left and right edges of the soil profile, a free drainage boundary condition at the bottom of the soil profile, and an atmospheric boundary condition at the top edge of the flow domain (Figure 4). In treatments SDI₁ and SDI₂, a time-variable flow boundary condition was defined at the nodes within a 6.75 cm-long segment at the top left of the flow domain, representing the water infiltration zone from the emitters. In treatments SSDI₁ and SSDI₂, a time-variable flow boundary condition was defined at the nodes within a 2 cm diameter hemisphere located 15 cm deep at the left lateral edge of the flow domain, representing the water emission zone in these treatments. Similarly, in treatment SSDI₃, a time-variable flow boundary condition was defined at the nodes within a 2 cm diameter hemisphere located 35 cm deep (25 cm deep in 2023) at the left lateral edge of the flow domain, representing the water emission zone in this treatment (Figure 4).

During the irrigation events, the irrigation flux was estimated by dividing the emitter discharge rate (1600 cm³ h⁻¹) by the surface wetted area. For treatments SDI₁ and SDI₂, this area was approximately 675 cm² (i.e., 50 cm × 13.5 cm, where the first value is the emitter spacing and the second is the width of the saturated zone through which water infiltrates from the emitter). This value was selected to ensure that irrigation water could infiltrate into the soil without producing positive surface pressure [34]. For treatments SSDI₁–SSDI₃, the irrigation flux was estimated by dividing the emitter discharge rate (1600 cm³ h⁻¹) by the flow area (50 cm × 6.28 cm, with the latter being the circumference length of the emitter pipe).

The initial soil moisture conditions differed between the two years studied. In 2022, due to the rainfall events that occurred prior to the planting date and continued during the first days after transplanting (Table 1), uniform moisture conditions equivalent to field capacity (matric potential = 200 cm H₂O) were defined across the entire flow domain. In 2023, due to the absence of rain in the months preceding planting, heterogeneous moisture conditions were established in the flow domain, starting from a matric potential of 1000 cm H₂O at the soil surface and decreasing to 100 cm H₂O at a depth of 100 cm (lower boundary of the flow domain).

In addition to the precipitation rate, the atmospheric boundary condition also uses as inputs daily potential transpiration (PT) and evaporation (PE) rates. Daily estimates of PT and PE were derived from daily crop evapotranspiration (ETc) values calculated using the FAO-56 approach [22] described in Section 2.2. Daily Kc values were derived from Equation (1), for which daily estimates of SS were obtained from the relationship between SS and cumulative growing degree days (GDDs) obtained in this study (Figure 5). GDDs were calculated using 10 °C and 30 °C as the lower and upper threshold temperatures, respectively, as proposed by Pathak and Stoddard [35] for this species. The partitioning of ETc into PT and PE assumed that PE gradually decreases with SS throughout the crop cycle to minimum values of 5% ETc for SS greater than or equal to 80%.

2.7. Water Footprint

According to the methodology proposed by Hoekstra et al. [5], the water footprint of a crop (WF) is the sum of three components:

WF = WF_green + WF_blue + WF_grey

(4)

and it is usually expressed in m³/t or L/kg.

The green component of the water footprint of a crop is calculated using the following expression:

$${WF}_{green}\left({\displaystyle \frac{m^3}{t}}\right)={\displaystyle \frac{{ET}_{green}\left(\frac{m^3}{ha}\right)}{Yield\left(\frac{t}{ha}\right)}}$$

(5)

where ET_green represents the fraction of rainfall that is used in the evapotranspiration process throughout the crop cycle.

The blue component of the water footprint of a crop is calculated using the following expression:

$${WF}_{blue}\left({\displaystyle \frac{m^3}{t}}\right)={\displaystyle \frac{{ET}_{blue}\left(\frac{m^3}{ha}\right)}{Yield\left(\frac{t}{ha}\right)}}$$

(6)

where ET_blue represents the fraction of irrigation water that is used in the evapotranspiration process throughout the crop cycle.

The grey component of the water footprint of a crop is calculated using the following expression:

$${WF}_{grey}\left({\displaystyle \frac{m^3}{t}}\right)={\displaystyle \frac{\frac{N_{excess}}{\left(C_{max}-C_{nat}\right)}\left(\frac{m^3}{ha}\right)}{Yield\left(\frac{t}{ha}\right)}}$$

(7)

where N_excess (kg ha⁻¹) represents the amount of nitrogen lost through leaching or runoff, C_max represents the maximum permissible concentration of N in the receiving water body, and C_nat is the existing concentration of N in the receiving water body before the contaminating activity occurs (considered negligible in this study). The value of C_max considered in this study is 11.29 mg N L⁻¹, corresponding to 50 mg $\mathit{N}{\mathit{O}}_{\mathbf{3}}^{-}$ L⁻¹ (EU Nitrates Directive, 91/676/EEC). Although nitrate in the soil solution is not the only substance that can contaminate groundwater and surface water sources, it is often considered the most critical. Determining the gray water footprint based on the most critical pollutant is considered sufficient as an overall indicator of water pollution [5]. Although the nitrate concentration measured in the soil solution may differ from what actually reaches freshwater sources due to adsorption processes, nitrate ion adsorption is low in soils poor in organic matter [36], such as the soil in this study.

The seasonal values of ET_green and ET_blue were obtained using the Hydrus 2D/3D program [26], as described previously. To differentiate between the fractions of ET that are green water and blue water, an initial simulation was performed without including irrigation inputs, which allowed for the determination of the ET_green component as the sum of the simulated actual evaporation and transpiration fluxes. Subsequent simulations included the volumes of irrigation water applied in each treatment, determining the ET_blue component as the difference between total ET and ET_green [5].

The simulated seasonal drainage flow was used, along with the nitrate concentrations measured in the soil solution (Section 2.3), to determine the amounts of nitrogen lost through deep percolation (N_excess).

2.8. Statistical Analyses

The effect of irrigation treatments on the various studied variables was determined using analysis of variance (ANOVA). The Shapiro–Wilks and Bartlett tests were used to check for normality and homoscedasticity of the data (p > 0.05). Mean separation was performed using Tukey’s test (p < 0.05). The coefficient of determination (R²) was utilized to assess the goodness of fit of the regression analyses conducted. The statistical analysis was conducted using the R statistical package [37].

3. Results

In 2022, the crop’s vegetative development was not affected by irrigation treatments (

Table 5. Fraction of the area shaded by the crop (%) during the 2022 and 2023 growth cycles. The values shown are means ± standard error. DAT: days after transplanting.

Growth Cycle 2022		Fraction of Shaded Area (%)
DAT	SDI1	SDI2	SSDI1	SSDI2	SSDI3	p-Value
7	0.2 ± 0.1	0.2 ± 0.1	0.2 ± 0.1	0.4 ± 0.2	0.4 ± 0.1	0.664
17	4.1 ± 2.4	2.1 ± 0.3	2.4 ± 0.2	5.6 ± 2.0	4.8 ± 1.4	0.313
24	11.1 ± 2.5	8.9 ± 1.6	12.2 ± 1.2	15.8 ± 2.9	13.5 ± 2.4	0.185
31	23.4 ± 0.7b	25.2 ± 2.4ab	36.1 ± 1.5a	29.1 ± 1.2ab	29.5 ± 1.3ab	0.006
38	44.3 ± 1.6	47.5 ± 2.5	49.4 ± 2.4	47.6 ± 2.0	46.3 ± 1.8	0.235
45	53.3 ± 4.1	55.4 ± 2.1	59.2 ± 4.4	56.5 ± 2.5	54.6 ± 2.5	0.628
52	68.6 ± 6.0	73.7 ± 2.7	75.2 ± 5.7	74.0 ± 4.4	70.6 ± 5.2	0.746
60	73.9 ± 4.2	82.4 ± 3.3	83.2 ± 4.9	78.6 ± 3.9	77.7 ± 5.2	0.635
Growth Cycle 2023
DAT	SDI1	SDI2	SSDI1	SSDI2	SSDI3	p-Value
11	0.7 ± 0.1	0.9 ± 0.1	0.6 ± 0.1	0.7 ± 0.1	0.8 ± 0.1	0.131
18	3.9 ± 0.3	4.6 ± 0.6	3.2 ± 0.4	4.2 ± 0.4	5.0 ± 0.5	0.079
25	15.7 ± 1.5	14.9 ± 1.5	14.9 ± 0.6	15.7 ± 1.5	18.4 ± 1.4	0.327
32	33.6 ± 2.4	34.6 ± 2.6	33.3 ± 2.1	36.2 ± 1.5	39.4 ± 1.3	0.220
39	49.3 ± 1.2b	49.4 ± 4.3b	49.9 ± 1.1ab	55.1 ± 1.3ab	59.0 ± 1.5a	0.023
46	58.7 ± 2.5b	60.2 ± 4.7ab	60.7 ± 3.0ab	65.1 ± 1.1ab	71.6 ± 2.1a	0.035
52	68.7 ± 3.0b	70.8 ± 6.3ab	76.1 ± 2.2ab	79.9 ± 2.6ab	79.8 ± 2.1a	0.041
63	77.6 ± 5.1	79.7 ± 6.1	88.7 ± 2.0	92.5 ± 1.6	91.4 ± 1.5	0.068

In each row, mean values followed by different letters indicate significant differences (p < 0.05).

). Significant differences in SS between treatments were observed on a single sampling date (DAT 31), with SDI₁ showing significantly lower SS values than SSDI₁ (Table 5). Two months after transplanting, all treatments showed similar SS values (p = 0.635), with crop coverage values close to 80% (Table 5), like those observed for this species under similar soil and climatic conditions [23]. In 2023, the irrigation treatments also showed similar vegetative development values (p = 0.068) at the end of the vegetative growth stage (DAT 63), although SSDI₃ showed higher values than SDI₁ during part of the vegetative growth phase.

The crop’s water status was significantly affected by irrigation treatments on three of the five sampling dates in 2022 (

Table 6. Leaf water potential (MPa) during the 2022 and 2023 growth cycles. The values shown are means ± standard error. DAT: days after transplanting; WSI: water stress integral.

Growth Cycle 2022		Leaf Water Potential (MPa)
DAT	SDI1	SDI2	SSDI1	SSDI2	SSDI3	p-Value
46	−0.67 ± −0.02	−0.66 ± −0.02	−0.71 ± −0.02	−0.66 ± −0.02	−0.62 ± −0.02	0.052
53	−0.76 ± −0.02a	−0.76 ± −0.02a	−0.82 ± −0.01ab	−0.74 ± −0.02a	−0.90 ± −0.05b	0.003
64	−0.85 ± −0.02b	−0.81 ± −0.01b	−0.79 ± −0.02b	−0.94 ± −0.02c	−0.71 ± −0.01a	<0.001
79	−0.76 ± −0.02	−0.80 ± −0.01	−0.75 ± −0.02	−0.77 ± −0.04	−0.71 ± −0.01	0.139
88	−1.06 ± −0.10ab	−0.93 ± −0.04ab	−0.98 ± −0.02ab	−1.14 ± −0.03b	−0.84 ± −0.02a	0.013
WSI	10.9 ± 0.8b	10.2 ± 0.3ab	10.3 ± 0.2ab	12.3 ± 0.3b	8.2 ± 0.5a	0.002
Growth Cycle 2023
DAT	SDI1	SDI2	SSDI1	SSDI2	SSDI3	p-Value
52	−0.49 ± −0.01	−0.51 ± −0.02	−0.51 ± −0.01	−0.50 ± −0.03	−0.47 ± −0.02	0.687
63	−0.89 ± −0.01b	−0.73 ± −0.03a	−0.78 ± −0.02ab	−0.78 ± −0.03a	−0.80 ± −0.03ab	0.008
73	−0.84 ± −0.02	−0.83 ± −0.04	−0.83 ± −0.04	−0.79 ± −0.01	−0.84 ± −0.04	0.846
84	−0.77 ± −0.01	−0.80 ± −0.04	−0.78 ± −0.02	−0.79 ± −0.04	−0.80 ± −0.02	0.916
94	−0.84 ± −0.02a	−0.83 ± −0.02a	−0.93 ± −0.07ab	−1.01 ± −0.02b	−0.90 ± −0.02ab	0.014
WSI	13.3 ± 0.3	11.7 ± 0.5	12.7 ± 0.4	12.9 ± 0.6	12.9 ± 0.5	0.233

In each row, mean values followed by different letters indicate significant differences (p < 0.05).

). Treatment SSDI₃ showed a higher water stress level at the end of the vegetative growth stage (DAT 53), while SSDI₂ tended to show lower LWP values during the fruit ripening stage (DAT 64 onwards). As a result, the treatments with the highest accumulated water stress level during the measurement period were SDI₁ and SSDI₂, with WSI values of 10.9 and 12.3 MPa·day, respectively. In 2023, significant differences between irrigation treatments were observed on only two of the five sampling dates, with SDI₁ showing the lowest LWP values at the end of the vegetative growth phase (DAT 63) and SSDI₂ the lowest LWP during the water deficit period (DAT 94) (Table 6). No significant differences in WSI were observed between treatments in 2023. Regarding leaf gas exchange, differences between treatments were observed on only two of the four sampling dates (

Table 7. Stomatal conductance (g_s) and leaf photosynthesis rate (A) measured during the 2022 growth cycle. The values shown are means ± standard error. DAT: days after transplanting.

DAT	SDI1	SDI2	SSDI1	SSDI2	SSDI3	p-Value
gs (mmol m−2 s−1)
60	751 ± 42	623 ± 12	687 ± 44	691 ± 31	623 ± 32	0.121
64	575 ± 23	506 ± 72	582 ± 32	533 ± 48	528 ± 18	0.698
79	647 ± 35ab	501 ± 28bc	742 ± 58a	625 ± 23ab	412 ± 17c	<0.001
88	861 ± 63	724 ± 92	893 ± 49	808 ± 75	712 ± 43	0.336
A (μmol m−2 s−1)
60	18.5 ± 0.6ab	19.2 ± 0.1ab	19.1 ± 0.5ab	20.0 ± 0.4a	17.5 ± 0.2b	0.008
64	18.6 ± 1.6	17.9 ± 0.5	17.2 ± 1.4	15.8 ± 1.6	17.9 ± 0.8	0.611
79	18.6 ± 0.6ab	18.5 ± 1.0ab	19.6 ± 0.3a	19.0 ± 0.8a	15.0 ± 1.3b	0.022
88	13.7 ± 0.4	14.6 ± 1.0	15.9 ± 1.3	16.1 ± 0.8	16.7 ± 0.5	0.142

In each row, mean values followed by different letters indicate significant differences (p < 0.05).

), with SSDI₃ showing the lowest g_s and A values on those dates.

Nitrate concentration in the soil solution (

Table 8. Nitrate concentration in soil solution during the 2022 and 2023 growth cycles. The values shown are means ± standard error. DAT: days after transplanting.

Growth Cycle	DAT	[NO3−] (mg L−1)
2022	46	195.8 ± 98.34
	88	745.9 ± 52.80
2023	52	352.7 ± 40.70
	73	247.5 ± 52.50
	84	362.5 ± 2.50

) presented high variability between the two sampling dates in 2022, with a resulting seasonal mean value of 470.8 mg L⁻¹. In 2023, there was greater stability in the nitrate concentrations measured in the soil solution, with a seasonal mean value of 320.9 mg L⁻¹. These mean values were used in the WF_gray estimates presented below.

Table 9. Yield and total soluble solids (TSS). The values shown are means ± standard error.

	Yield (t ha−1)		TSS (°Brix)
Treatment	2022	2023	2022	2023
SDI1	128 ± 13	127 ± 11	5.5 ± 0.1	6.1 ± 0.2
SDI2	123 ± 8	135 ± 13	5.5 ± 0.2	6.3 ± 0.1
SSDI1	149 ± 6	150 ± 17	5.2 ± 0.2	5.7 ± 0.2
SSDI2	120 ± 4	157 ± 12	5.2 ± 0.2	5.9 ± 0.1
SSDI3	124 ± 3	160 ± 3	5.2 ± 0.5	5.7 ± 0.1
p-value	0.108	0.354	0.418	0.152

shows the yield and fruit quality results obtained in both cultivation cycles. In 2022, treatment SSDI₁ had the highest average yield (149 t/ha), while treatments SSDI₁–SSDI₃ had the highest average yields (150–160 t/ha) in 2023. Due to the high variability observed between replications, the observed differences were not statistically significant. The concentration of total soluble solids (TSS) expressed as degrees Brix was very similar across all treatments, ranging from 5.2 to 5.5 in 2022 and from 5.7 to 6.3 in 2023.

Table 10. Cumulative values of the components of the soil water balance throughout the entire growth cycle. P: precipitation; I: irrigation; PE: potential evaporation; PT: potential transpiration; AE: actual evaporation; AT: actual transpiration; DD: deep drainage.

Growth Cycle 2022
T	P (mm)	I (mm)	PE (mm)	PT (mm)	AE (mm)	AT (mm)	DD (mm)
SDI1	57.4	570.0	137.0	497.0	114.4	490.8	35.6
SDI2	57.4	561.5	137.0	497.0	111.5	490.8	33.0
SSDI1	57.4	517.3	137.0	497.0	108.4	462.2	34.8
SSDI2	57.4	495.3	137.0	497.0	108.7	449.1	35.2
SSDI3	57.4	565.9	137.0	497.0	84.4	391.4	126.6
Growth Cycle 2023
T	P (mm)	I (mm)	PE (mm)	PT (mm)	AE (mm)	AT (mm)	DD (mm)
SDI1	98.4	494.1	122.5	408.8	103.0	406.9	44.3
SDI2	98.4	465.2	122.5	408.8	101.8	401.4	35.0
SSDI1	98.4	412.9	122.5	408.8	100.7	379.7	22.3
SSDI2	98.4	385.3	122.5	408.8	100.6	357.7	21.6
SSDI3	98.4	421.8	122.5	408.8	91.6	365.6	40.2

shows the simulated seasonal values of the soil water balance components. The amount of water lost via soil evaporation in SDI treatments (SDI₁, SDI₂) represented 82–84% of potential soil evaporation in both cultivation cycles. SSDI treatments with shallowly buried driplines (SSDI₁ and SSDI₂) showed soil evaporation losses very similar to those of SDI₁ and SDI₂, representing 79% (2022) and 82% (2023) of potential soil evaporation. Treatment SSDI₃, whose driplines were buried deeper, showed a significantly lower amount of evaporated water from the soil surface than the other treatments, representing 61% (2022) and 75% (2023) of potential soil evaporation.

Actual transpiration represented, in 2022, 99% (SDI₁), 99% (SDI₂), 93% (SSDI₁), 90% (SSDI₂), and 79% (SSDI₃) of potential transpiration, while in 2023, it represented 99% (SDI₁), 98% (SDI₂), 93% (SSDI₁), 87% (SSDI₂), and 89% (SSDI₃) of potential transpiration. In 2022, although the amount of irrigation applied to SSDI₃ was very similar to that applied to SDI₁ and SDI₂ (Table 10), actual transpiration in SSDI₃ was 20% lower than in these treatments, with SSDI₃ showing a significantly higher amount of water lost via deep drainage. In 2023, with a slightly shallower installation depth of the driplines, the amount of irrigation applied in SSDI₃ was 15% lower than that applied in SDI₁, while the amount of water transpired in SSDI₃ was 10% lower than in SDI₁. Unlike what was observed in 2022, the water losses via drainage in SSDI₃ in 2023 were lower than in 2022 and similar to those observed in SDI₁, confirming that a dripline depth of 25 cm is a more suitable option than a depth of 35 cm for the soil conditions and crop species evaluated in this study.

Table 11. Green, blue, gray, and total water footprints (WF) calculated for all irrigation treatments and the 2022 and 2023 growth cycles.

Growth Cycle 2022
Treatment	WFgreen (L kg−1)	WFblue (L kg−1)	WFgrey (L kg−1)	WFtotal (L kg−1)
SDI1	8.1	39.1	26.1	73.3
SDI2	8.4	40.4	25.1	73.9
SSDI1	7.0	31.3	22.0	60.3
SSDI2	8.7	37.9	27.7	74.3
SSDI3	8.4	29.9	95.6	133.8
Growth Cycle 2023
Treatment	WFgreen (L kg−1)	WFblue (L kg−1)	WFgrey (L kg−1)	WFtotal (L kg−1)
SDI1	8.5	31.7	22.4	62.6
SDI2	8.0	29.3	16.7	54.0
SSDI1	7.2	24.9	9.6	41.6
SSDI2	6.8	22.3	8.8	38.0
SSDI3	6.7	21.9	16.2	44.8

shows the water footprint obtained for the five treatments during the two cultivation cycles evaluated. The green water footprint (WF_green) was very similar across all treatments, ranging from 7.0 (SSDI₁) to 8.7 (SSDI₂) L kg⁻¹ in 2022 and from 6.7 to 8.5 L kg⁻¹ in 2023 (Table 11). Larger differences were observed in the blue water footprint (WF_blue) values, both between treatments and between cultivation cycles. In 2022, WF_blue ranged from 29.9 (SSDI₃) to 40.4 (SDI₂) L kg⁻¹, while in 2023, it ranged from 21.9 (SSDI₃) to 31.7 (SDI₁) L kg⁻¹. The gray water footprint (WF_gray) values also varied greatly between the two cultivation cycles. In 2022, WF_gray was similar among treatments SDI₁–SSDI₂, ranging from 22.0 (SSDI₁) to 27.7 (SSDI₂) L kg⁻¹. However, treatment SSDI₃ showed a substantially higher WF_gray (95.6 L kg⁻¹) than the other treatments due to the reported high deep drainage depths observed in SSDI₃ (Table 10). In 2023, all treatments showed lower WF_gray than in 2022, ranging from 8.8 (SSDI₂) to 22.4 (SDI₁) L kg⁻¹. The lower nitrate concentration observed in the soil solution in 2023, along with the lower amount of percolated water, especially in SSDI treatments, would explain these results.

In terms of total water footprint (WF_total), SSDI₁ showed the lowest value in 2022 (60.3 L kg⁻¹), which was 18% lower than the values observed in SDI₁, SDI₂, and SSDI₂, and 55% lower than that of SSDI₃. In 2023, WF_total values were lower than those obtained in 2022, with treatment SSDI₂ showing the lowest WF_total (38 L kg⁻¹) and SDI₁ the highest WF_total (62.6 L kg⁻¹). In percentage terms and using aggregated data from the two years of study (Figure 6), the WF_green represented approximately 12–14% of WF_total in SDI₁–SSDI₂, and around 8.5% in SSDI₃. WF_blue represented just over 50% of WF_total in SDI₁–SSDI₂, representing around 30% of WF_total in SSDI₃. WF_gray accounted for approximately 30–35% of WF_total in SDI₁–SSDI₂, rising to approximately 60% in treatment SSDI₃.

4. Discussion

The evaluated irrigation treatments did not significantly impact the crop’s vegetative development (Table 5). This aligns with previous studies, such as Patane et al. [19], who observed that significant reductions in total crop biomass only occurred when ET reductions exceeded 50% compared to full irrigation. In this study, the lowest ET values were observed in treatment SSDI₃ (2022) and treatments SSDI₂ and SSDI₃ (2023), with ET reductions of approximately 20% in 2022 and 10% in 2023 compared to treatment SDI₁ (Table 10). These reductions were not severe enough to affect the vegetative development of the crop, consistent with the findings of Patane et al. [19].

The close relationship between SS and GDDs over the two years of study (Figure 5) enabled the determination of daily water needs, essential input data for the Hydrus 2D/3D model to compute soil water balance components. This relationship is also valuable for developing precision irrigation tools, allowing the prediction of the daily crop coefficient (Equation (1)) based solely on thermal data. Additionally, it helps identify deviations in the crop’s vegetative growth pattern through simple, periodic SS measurements during the vegetative growth phase.

The experimental results indicate that SSDI is a promising alternative to SDI, potentially enhancing the crop’s environmental sustainability. Treatment SSDI₁, with a 10% irrigation reduction in 2022 and 15% in 2023 compared to SDI₁, showed no water stress symptoms (Table 6 and Table 7) nor reductions in fruit yield and quality (Table 9). These findings suggest that the Kc values currently available in the literature—often based on single crop coefficients—should be updated to better reflect local cultivation conditions. It would be necessary to conduct lysimeter-based studies to determine precise dual crop coefficients [22] adjusted to local conditions, especially valuable for SSDI irrigation systems. Despite the greater water deficit during the fruit ripening phase and the tendency for lower LWP values during this period (Table 6), treatment SSDI₂ performed similarly to the other treatments and saved 13% (2022) and 22% (2023) of irrigation water compared to SDI₁. Previous studies by Patane et al. [19] and Del Amor and Del Amor [38] also indicated that severe ET reductions were required to affect tomato yield significantly. In this study, numerical simulations with the Hydrus model (Table 10) showed that the fraction of water lost by evapotranspiration in SSDI₂ was 8% (2022) and 10% (2023) lower than that of SDI₁, well below the productivity loss thresholds observed in other studies [22].

The depth at which driplines are buried in SSDI treatments is crucial for water application efficiency. According to Lamm [20], common depths for processing tomato cultivation range from 0.15–0.35 m. Numerical simulations confirmed that a depth of 35 cm (SSDI₃ in 2022) led to higher water drainage losses than shallower depths (Table 10), resulting in lower LWP and leaf photosynthesis values on some 2022 sampling dates (Table 6 and Table 7), although these did not translate into productivity losses (Table 9). Reducing the dripline depth in treatment SSDI₃ to 25 cm (2023) brought drainage losses closer to those of shallower installations (Table 10).

A significant interannual influence was observed on the crop’s water footprint (Table 11). WF_total ranged from 60.3 L kg⁻¹ (SSDI₁) to 133.8 L kg⁻¹ (SSDI₃) in 2022, and from 38.0 L kg⁻¹ (SSDI₂) to 62.6 L kg⁻¹ (SDI₁) in 2023, driven mainly by differences in evaporative demand (ET₀: 677.8 mm in 2022 vs. 579.5 mm in 2023, Table 1). Analyzing the water footprint across seasons with significant variations in evaporative demand was crucial in this work, as this has proven to be a key factor when using the water footprint as a comparative parameter between different production areas or seasons. This variability is important to consider when estimating water footprints for agricultural production, an increasingly frequent practice due to the growing interest in including WF as a sustainability criterion in certifications and eco-labels of agri-food products, which can positively influence consumer purchasing decisions. The irrigation system used also significantly affects WF estimates. In both growing cycles, SSDI had significantly lower WF_total values than the control treatment (SDI₁). In 2022, SSDI₁ had a WF_total approximately 20% lower than SDI₁, and in 2023, WF_total values for SSDI₁ and SSDI₂ were 35% lower than for SDI₁.

Previous studies estimating the WF of tomato cultivation in different Mediterranean regions show considerable disparity in WF values, likely due to differences in edaphoclimatic conditions and the uncertainty in WF_gray estimates. Chico et al. [12] determined the WF of tomato production across Spain, with average WF_total values of 236 L kg⁻¹ (5 L kg⁻¹ WF_green, 92 L kg⁻¹ WF_blue, and 139 L kg⁻¹ WF_gray). Except for WF_green, these values are significantly higher than those observed in this study (Table 11). Chapagain and Orr [39] reported average WF_total values of 81.3 L kg⁻¹ for Spain and 70.1 L kg⁻¹ for Andalusia, closer to (though still higher than) those in this study. Aldaya and Hoekstra [3] found significantly higher WF_total values (114 L kg⁻¹) for processing tomatoes in Italy, mainly due to higher WF_blue values (60 L kg⁻¹) compared to this study. Conversely, the WF_total of processing tomatoes in central Greece (61 L kg⁻¹) is similar to some values observed in this study (e.g., SSDI₁ in 2022: 60.3 L kg⁻¹; SDI₁ in 2023: 62.6 L kg⁻¹), with similar green, blue, and grey WF values.

Given the excellent performance of processing tomatoes under SSDI, including its combination with deficit irrigation, future research should focus on performing a comprehensive economic analysis comparing SSDI and SDI irrigation systems. This analysis should consider installation costs, maintenance, water savings, and crop yield to provide a clearer picture of the financial implications and benefits for farmers. Additionally, investigating the factors influencing farmer adoption of SSDI is essential for developing policy recommendations that encourage its uptake. This includes understanding barriers to adoption and providing necessary training and resources to farmers. These research directions will help to better understand the economic viability and practical challenges of SSDI, facilitating its broader implementation and promoting sustainable agricultural practices.

5. Conclusions

The results obtained in this study demonstrate that subsurface drip irrigation (SSDI) makes more efficient use of water resources than surface drip irrigation (SDI). With a dripline depth of 15 cm and a water application rate 10–15% lower than that applied to the control treatment (surface drip irrigation with no water limitation, SDI₁), the SSDI₁ treatment did not show water stress symptoms or yield reductions, reducing the total water footprint by 20–35% compared to SDI₁. These findings suggest the need for lysimeter-based studies or other approaches to determine crop coefficients (ideally dual Kc) for processing tomatoes grown under SSDI and for the local conditions in Andalusia. SSDI combined with a mild-to-moderate water deficit (15–20%) during the fruit ripening phase did not significantly reduce yield and showed similar (2022) or markedly lower (2023) water footprint values compared to SDI₁. When the dripline depth was 35 cm, SSDI₃ exhibited a much higher water footprint than the other SDI and SSDI treatments due to a high WF_grey caused by significant drainage losses. However, when the dripline depth was reduced to 25 cm in the second year of the study, the performance of SSDI₃ was very similar to that of SSDI₁. Inter-seasonal differences in atmospheric water demand caused significant variations in total water footprint (lower WF with reduced evaporative demand), highlighting its importance as a key factor when using water footprint as a comparative indicator between different production areas or seasons. Methodologically, this study also highlighted the importance of including all components of the water footprint in the calculation of the total water footprint, as WF_grey (often overlooked) can be decisive in scenarios where deep percolation losses are high. This study, due to the lack of recent precedents in Andalusia, will serve as a baseline for local producers and for decision-making in irrigation management based on modern precision irrigation techniques.