Modeling Approaches for Determining Dripline Depth and Irrigation Frequency of Subsurface Drip Irrigated Rice on Different Soil Textures

Water saving techniques such as drip irrigation are important for rice (Oriza sativa L.) production in some areas. Subsurface drip irrigation (SDI) is a promising alternative for intensive cropping since surface drip irrigation (DI) requires a higher degree of labor to allow the use of machinery. However, the semi-aquatic nature of rice plants and their shallow root system could pose some limitations. A major design issue when using SDI is to select the dripline depth to create appropriate root wetting patterns as well as to reduce water losses by deep drainage and evaporation. Soil texture can greatly affect soil water dynamics and, consequently, optimal dripline depth and irrigation frequency needs. Since water balance components as deep percolation are difficult to estimate under field conditions, soil water models as HYDRUS-2D can be used for this purpose. In the present study, we performed a field experiment using SDI for rice production with Onice variety. Simulations using HYDRUS-2D software successfully validated soil water distribution and, therefore, were used to predict soil water contents, deep drainage, and plant water extraction for two different dripline depths, three soil textures, and three irrigation frequencies. Results of the simulations show that dripline depth of 0.15 m combined with one or two daily irrigation events maximized water extraction and reduced percolation. Moreover, simulations with HYDRUS-2D could be useful to determine the most appropriate location of soil water probes to efficiently manage the SDI in rice.


Introduction
Rice is a main international crops and is a staple food in many countries. Significantly, 90% of worldwide rice production takes place in Asia. In 2018, 67% of the world's total rice production came from China, India, Indonesia, and Bangladesh [1]. In Europe, Spain is the second largest producer after Italy, with a total surface area of 105,422 ha in 2019 [2].
About 75% of rice is produced under continuous flooded irrigation and uses 34 to 43% of the world's irrigation water [3]. The total average water use in continuously flooded irrigation can greatly vary depending on the soil type, rice culture, and water management practices, ranging from 6750 m 3 ha −1 to 44,500 m 3 ha −1 [4].
Water scarcity represents a growing concern in vast producing areas of the world [5]. For this reason, the use of water saving irrigation techniques becomes critical. Alternate wetting and drying irrigation management technique and aerobic rice production systems have shown its effectiveness three-dimensional (3D) approach. Generally, when emitters are close to each other along the dripline an infinite line source is considered and a 2D model is applied [37][38][39].
Reyes-Esteves and Slack [34] used HYDRUS-2D to simulate spatial and temporal distributions of soil water contents, root water uptake of an SDI system commonly used for alfalfa production in California. In this work, the depth at which driplines must be placed in alfalfa was optimized considering different scenarios of irrigation patterns and soil textures. At the knowledge of the authors there are no previous applications of HYDRUS-2D on rice production using subsurface drip irrigation.
The general objectives of this work were to determine the depth at which driplines should be installed, to establish the most effective irrigation frequency for rice production using SDI in soils of different textures, and to estimate major components of the field water balance. Field experiments and simulations with HYDRUS-2D were carried out for these purposes.

Field Experiment
Field tests were carried out in Pals (Baix Ter, Girona, Spain) during 2019 growing season on an experimental plot of 0.39 ha (41 • 58 43" N, 3 • 08 21" E). The soil texture of the plot was spatially heterogeneous; sandy-loam textural class prevailed in 85% of the surface area of the experimental field and loamy one in the other 15% (Table 1). In addition, another soil of silty-clay texture was characterized (Table 1) since it was also representative of the Baix Ter production area, with the further objective of extrapolating the results obtained in the field test with the HYDRUS-2D simulation results. The SDI system was installed with driplines buried at 0.15 m and spaced 0.66 m. Thin-walled driplines (0.38 mm), 16.20 mm internal diameter with integrated pressure-compensating emitters of 1 L h −1 (Dripnet PC TM AS, Netafim, Israel), spaced 0.3 m apart, were used.
The plot was sowed on 28 May 2019 using the variety Onice on dry soil conditions with a seed drill at 135 kg ha −1 and spacing the seeding lines 0.19 m. The general irrigation criteria was based on maintaining −10 kPa of soil water potential at 0.10 m depth and at a horizontal distance of 0.08 m from the emitter (Figure 1). This matric potential corresponds to the field capacity for aerobic rice [40].
Temperature, rainfall, atmospheric pressure, air humidity, global radiation, and wind speed were measured with an automated meteorological station, located approximately 2 km from the experimental plot. The reference evapotranspiration (ET o ) was determined using the FAO-56 method [41], which requires solar radiation, maximum, minimum, and mean air temperatures, air humidity, and wind speed data measured with the automated meteorological station.
Irrigation water applied to the plot was continuously measured using a volumetric meter. Soil water contents were continuously determined from day 34 after sowing (DAS), at 0.10 and 0.25 m depths and in two horizontal positions, one close to the emitter, at 0.08 m from the dripline, and the other at the intermediate region between the two driplines ( Figure 1). Volumetric soil water contents were measured using a sensor based on frequency domain reflectometry, with an accuracy of ±3% (SM 150-T, Delta-T Devices, Cambridge, UK). Sensors were located in part of the plot with a sandy-loam texture. Depth of the water table was monitored with weekly readings in 3 piezometers, 2 m long, located around the plot, which allowed to establish adequate boundary conditions in simulations, as explained in Section 2.2.
Phenological stages of the crop were monitored throughout the campaign. Yield was obtained from sampling carried out on 18 October 2019, at the stage of physiological maturity. Six samples were taken on a reference surface of 0.5 m × 0.5 m for each of the two soil textures present in the plot (sandy-loam and loam, Table 1). After drying the samples, standardized production was obtained at a grain humidity of 14%. Average productions in the surfaces corresponding to each textural class were determined, carrying out an analysis of variance to find out the existence of significant differences depending on soil texture. piezometers, 2 m long, located around the plot, which allowed to establish adequate boundary conditions in simulations, as explained in Section 2.2.
Phenological stages of the crop were monitored throughout the campaign. Yield was obtained from sampling carried out on 18th October 2019, at the stage of physiological maturity. Six samples were taken on a reference surface of 0.5 m × 0.5 m for each of the two soil textures present in the plot (sandy-loam and loam, Table 1). After drying the samples, standardized production was obtained at a grain humidity of 14%. Average productions in the surfaces corresponding to each textural class were determined, carrying out an analysis of variance to find out the existence of significant differences depending on soil texture.
Water productivity (WP), which characterizes yield (Y) per unit of water used, was calculated in two ways: being Y the rice grain yield at 14% of humidity (kg ha −1 ); Irr the irrigation water (m 3 ha −1 ); and R the precipitation during the cultivation period (m 3 ha −1 ).  Water productivity (WP), which characterizes yield (Y) per unit of water used, was calculated in two ways: being Y the rice grain yield at 14% of humidity (kg ha −1 ); Irr the irrigation water (m 3 ha −1 ); and R the precipitation during the cultivation period (m 3 ha −1 ).

Simulations of Soil Water Dynamics
In order to model soil water dynamics, HYDRUS-2D [42] was selected. Due to the short distance between two consecutive emitters along the dripline (0.3 m), soil water movement can be assumed as two-dimensional. This well-known software package solves numerically Richards equation: Water 2020, 12, 1724

of 15
where θ is the volumetric soil water content (m 3 m −3 ), t is the time (day), x and z the spatial horizontal and vertical coordinates (m), respectively, h the soil water pressure (m), K(h) the hydraulic conductivity function (m day −1 ), and S(h) the sink term in the equation (m 3 m −3 day −1 ), to consider plant water extraction. Soil water retention and hydraulic conductivity functions were determined following the van Genuchten-Mualem model [43]. The parameters of the equations are shown in Table 2 and were estimated with Rosetta [44] from the percentages of sand, silt, and clay. Bulk density and water content were at −33 and −1500 kPa (Table 1). A flow domain 1 m deep and 0.33 m wide was considered, going from the vertical line where the dripline was located to the intermediate position between two driplines. Spatial discretization was performed with triangular elements with a dimension close to 0.001 m in the area closest to the emitter up to approximately 0.020 m in the area furthest away. In all simulations, initial condition in the whole domain was fixed as a matric potential of −33 kPa, corresponding to the initial conditions in the field test in the sandy-loam textured soil. In order to compare the results with the ones of the rest of simulated scenarios, the same initial condition was used. The boundary conditions were absence of flow on the two sides of the domain, unit vertical hydraulic gradient simulating free drainage from the bottom of the domain since, according to the piezometer readings, the depth of the aquifer was always greater than 2 m, and an atmospheric boundary condition at the soil surface. In the elements representing the emitter, a variable flow condition was considered over time, in which the flow-rate during irrigation periods was obtained according to: being F, the water-flow applied to the variable condition in the emitter zone (m day −1 ); q emitter discharge (m 3 day −1 ); d emitter spacing (m); and r external radius of the dripline (m). In addition to the simulated field tests, different scenarios collected in the A-V cases listed in Table 4 were considered: two dripline depths (0.15 and 0.25 m) and three different soil textures (silty-clay, loam, and sandy-loam) that were representative of Baix Ter rice production area. It should be noted that the depths considered for the driplines are compatible with the cultural tasks foreseen.
Simulated irrigation practices were: (1) irrigation practiced during 2019 campaign (cases A-C, Table 4); (2) irrigation based strictly on crop evapotranspiration applying three different irrigation frequencies: two irrigations per day, one daily irrigation, and one irrigation every four days (Cases D-V, Table 4). A conceptual figure showing input and output parameters and followed procedure for the set of simulations carried out is depicted in Figure 2.
Climatic data from sowing (18 May 2019) to harvest (18 October 2019) was considered. Crop evapotranspiration (ET c ) was determined to consider water extraction in simulations according to the equation: being K c the crop coefficient, determined for the reference region following the suggested values by Allen et al. [41]. evapotranspiration (ETc) was determined to consider water extraction in simulations according to the equation: being Kc the crop coefficient, determined for the reference region following the suggested values by Allen et al. [41]. HYDRUS-2D code requires evaporation and transpiration as two independent input terms, since evaporation is treated as a surface boundary condition and transpiration as the sink term in the Richards equation. For a crop that does not cover the entire surface, which happened from sowing time up to 70 days later, the potential evapotranspiration is divided between potential evaporation (Ep) and potential transpiration (Tp). The partition of these two terms was carried out using the leaf area index (LAI), following the procedure described in the SWATRE model and using the equations [45]: Global solar radiation coefficient in rice (Kgr) was set to 0.3 according to Li et al. [17] and Phogat et al. [46]. LAI values were determined based on the phenological stages observed in the same plot and with the same variety of rice during the 2018 irrigation season. LAI values were obtained every two weeks using PocketLAI [47]. During the simulated periods in which there was no crop and potential evaporation was considered to be equal to ETo.
Water extracted by the crop was determined according to the equations proposed by Feddes [48] and with the optimized response parameters for rice crop obtained by Singh et al. [49] and applied, among others, by Li et al. [50]. It was considered that water extraction was located at the depth where the rice root system develops, being more intense near the soil surface and decreasing towards the maximum rooting depth, which was 0.25 m, according to field observations throughout the irrigation season. HYDRUS-2D code requires evaporation and transpiration as two independent input terms, since evaporation is treated as a surface boundary condition and transpiration as the sink term in the Richards equation. For a crop that does not cover the entire surface, which happened from sowing time up to 70 days later, the potential evapotranspiration is divided between potential evaporation (E p ) and potential transpiration (T p ). The partition of these two terms was carried out using the leaf area index (LAI), following the procedure described in the SWATRE model and using the equations [45]: Global solar radiation coefficient in rice (K gr ) was set to 0.3 according to Li et al. [17] and Phogat et al. [46]. LAI values were determined based on the phenological stages observed in the same plot and with the same variety of rice during the 2018 irrigation season. LAI values were obtained every two weeks using PocketLAI [47]. During the simulated periods in which there was no crop and potential evaporation was considered to be equal to ET o .
Water extracted by the crop was determined according to the equations proposed by Feddes [48] and with the optimized response parameters for rice crop obtained by Singh et al. [49] and applied, among others, by Li et al. [50]. It was considered that water extraction was located at the depth where the rice root system develops, being more intense near the soil surface and decreasing towards the maximum rooting depth, which was 0.25 m, according to field observations throughout the irrigation season.

Validation of the Simulation Model
HYDRUS-2D model was validated by comparing simulated soil water contents with soil water contents measured in the experimental plot at the positions described in Section 2.1. The validation period included the whole irrigation campaign and after the harvest, from 1 July 2019, when sensors were installed, until 22 January 2020. The root of the mean square error (RMSE), normalized RMSE (nRMSE), the coefficient of determination (R 2 ) and Nash-Sutcliffe model efficiency coefficient (E) were calculated to quantify the goodness of model predictions as follows: where P i are simulated values, O i are observed values, O is the mean of observed values, P is the mean of predicted values, and n is the total number of observations. The RMSE indicates the difference between observed and simulated soil water contents. It has the advantage of having the same units as measured value and therefore it can be compared with it. A lower value of RMSE indicates a good prediction of the model, meaning a value of zero is a perfect fit. The nRMSE relativizes the root mean square error value, dividing it by the mean of the observed values. When nRMSE is less than 10%, the prediction of the model is excellent, whereas the prediction is good if its value ranges between 10-20%. Further, the prediction is fair if it is greater than 20% and lower than 30%, and poor when it is greater than 30% [51,52]. E means it has been extensively used to evaluate the performance of hydraulic models, ranging its values from −∞ to 1, with greater values indicating better agreement. When E equals zero the observed mean O is as good as a predictor as the model. When E is negative it indicates that O is a better predictor than the model [53,54].

Water Balance Components
Simulation of soil water dynamics in the experimental plot and in the considered scenarios allowed the assessment of the water balance throughout the cultivation period (from the day of planting to harvesting, 144 days in total) considering water inlets (rain (R) and irrigation (Irr)) as well as outlets (drainage below 1 m and crop evapotranspiration). In addition, variation of soil water storage in the entire soil profile, error in the water balance of the model and the water that drained below 0.3 m depth (which was considered not usable by the crop, due to the depth of the rice rooting system) were computed. The relationship between cumulative evapotranspiration computed by the model (ET sim ) and water inlets in the soil profile was also determined. (ET sim / (Irr + R))·100 shows the percentage of water evapotranspirated in relation to contributions of water from irrigation and rain (Irr + R), being desirable to maximize its value to increase water use efficiency.

Irrigation Campaign in the Experimental Plot
Suitability of the irrigation pattern followed during the growing season is analyzed in this section. The amount of water supplied by irrigation during the first three weeks was higher than evapotranspiration demand (Figure 2). The reason for this was the need to promote seed germination after sowing. From day 15 to 70 DAS (from 12 June to 6 August 2019), the plot was watered in order to maintain the soil water content in sensor 2 (Figure 1) at 0.257 m 3 m −3 , corresponding to a matric potential of −10 kPa, resulting in cumulative applied irrigation water was lower than cumulative crop water demand (ET cc ) (Figure 3). To correct this problem, the irrigation schedule was modified from day 70 to 123 and during this period a slightly higher irrigation water amount than that of the extractions estimated from the ET c was provided. Consequently, from day 84 the accumulated irrigation surpassed the cumulative ET c again ( Figure 3). From day 123 until day 140, irrigation was stopped as the crop reached physiological maturity and the goal was to reduce the grain humidity in order to harvest as soon as possible.
section. The amount of water supplied by irrigation during the first three weeks was higher than evapotranspiration demand (Figure 2). The reason for this was the need to promote seed germination after sowing. From day 15 to 70 DAS (from June 12th to August 6th 2019), the plot was watered in order to maintain the soil water content in sensor 2 (Figure 1) at 0.257 m 3 m −3 , corresponding to a matric potential of −10 kPa, resulting in cumulative applied irrigation water was lower than cumulative crop water demand (ETcc) (Figure 3). To correct this problem, the irrigation schedule was modified from day 70 to 123 and during this period a slightly higher irrigation water amount than that of the extractions estimated from the ETc was provided. Consequently, from day 84 the accumulated irrigation surpassed the cumulative ETc again (Figure 3). From day 123 until day 140, irrigation was stopped as the crop reached physiological maturity and the goal was to reduce the grain humidity in order to harvest as soon as possible.   Figure 4. The matric potential at 0.08 m from the dripline and at 0.10 m depth was found above −10 kPa (corresponding to a volumetric content of 0.257 m 3 m −3 ) during the entire growing period (from 34 to 140 days DAS). Therefore, it should not suppose a reduction of potential yield according to Bouman et al. (2005) [40], as this water potential corresponds to field capacity for aerobic rice. At the same distance from the dripline but at a depth of 0.25 m, the matric potential was below −20 kPa (corresponding to 0.21 m 3 m −3 ) practically throughout the entire crop cycle. With this value of matric potential, notable reductions in production were expected according to Cabangon et al. (2003) [55].
In the furthest measurement points from the dripline (positions 3 and 4), the matric potential was also below −20 kPa between 34 and 70 DAS, corresponding to the period where irrigation water supply was lesser than ET c , meaning that the crop might be seriously stressed in the area located between driplines.

Validation of the Model
The visual comparison of the simulated and observed soil water contents in the experimental plot showed a similar trend (Figure 4), which illustrates a good performance of the model. Root mean square error (RMSE) values ranging from 0.028 m 3 m −3 and 0.064 m 3 m −3 , nRMSE ranged from 9.4% to 25.3%, determination coefficients (R 2 ) ranging from 0.38 and 0.54 and E from −0.589 to 0.353 (Table 3). Only in the positon (0.08, 0.10) the value of E was negative, meaning that the observed mean is a better predictor than the model. Nevertheless, the nRMSE in this case was less than 10% which indicates and excellent performance of the model. The greatest differences were found from 30 to 100 days after sowing in the area furthest from the dripline, where the observed soil water contents were systematically higher than that simulated with HYDRUS-2D.  Soil water contents measured at different positions from the dripline (Figure 1) are shown in Figure 4. The matric potential at 0.08 m from the dripline and at 0.10 m depth was found above −10 kPa (corresponding to a volumetric content of 0.257 m 3 m −3 ) during the entire growing period (from 34 to 140 days DAS). Therefore, it should not suppose a reduction of potential yield according to Bouman et al. (2005) [40], as this water potential corresponds to field capacity for aerobic rice. At the same distance from the dripline but at a depth of 0.25 m, the matric potential was below −20 kPa (corresponding to 0.21 m 3 m −3 ) practically throughout the entire crop cycle. With this value of matric potential, notable reductions in production were expected according to Cabangon et al. (2003) [55].
In the furthest measurement points from the dripline (positions 3 and 4), the matric potential was also below −20 kPa between 34 and 70 DAS, corresponding to the period where irrigation water supply was lesser than ETc, meaning that the crop might be seriously stressed in the area located between driplines.

Validation of the Model
The visual comparison of the simulated and observed soil water contents in the experimental plot showed a similar trend (Figure 4), which illustrates a good performance of the model. Root mean square error (RMSE) values ranging from 0.028 m 3 m −3 and 0.064 m 3 m −3 , nRMSE ranged from 9.4% to 25.3%, determination coefficients (R 2 ) ranging from 0.38 and 0.54 and E from −0.589 to 0.353 (Table  3). Only in the positon (0.08, 0.10) the value of E was negative, meaning that the observed mean is a better predictor than the model. Nevertheless, the nRMSE in this case was less than 10% which indicates and excellent performance of the model. The greatest differences were found from 30 to 100 days after sowing in the area furthest from the dripline, where the observed soil water contents were systematically higher than that simulated with HYDRUS-2D.
The statistics of comparison are similar to those obtained in previous works in which the observed and simulated soil water contents were compared in bare soil [37], in an apple orchard [16], and in laboratory scale experiments [56]. The differences between simulated and measured soil water contents can be explained by different causes, being some of them: (1) the soil variability, (2) the different representative volume of soil explored by the soil water sensor (which has a rod length of  The statistics of comparison are similar to those obtained in previous works in which the observed and simulated soil water contents were compared in bare soil [37], in an apple orchard [16], and in laboratory scale experiments [56]. The differences between simulated and measured soil water contents can be explained by different causes, being some of them: (1) the soil variability, (2) the different representative volume of soil explored by the soil water sensor (which has a rod length of 51 mm) and the nodal points where HYDRUS-2D provided the water content, and (3) the accuracy of the measurement techniques.
A calibration of the model with an independent dataset of measurements in the same soil can be used to reduce the differences between measured and simulated soil water contents [18,57]. Nevertheless, calibration requires an independent dataset of soil water measurements obtained in similar conditions which is not always feasible. Table 4 shows the variables that determine water balance terms, simulated cumulative evapotranspiration with HYDRUS-2D (ET sim ), percolation below 0.3 m, error and some derived relationships for each one of the analyzed cases, as described in Section 2.4. Table 4. Simulated water balance components (mm) considering 1 m depth of soil profile, water balance error (mm), relationship between simulated cumulative crop evapotranspiration (ET sim ), and water inlets, and drainage below the maximum depth of roots (0.3 m). Irr: Irrigation, R: Rain, ET sim : Simulated Evapotranspiration, P 1 m : Percolation below 1 m depth; SS: Soil storage change (negative when the content increases), BE: Total water balance error, ET sim /(Irr + R): Relationship in percentage between simulated crop evapotranspiration and soil water inlet (Irr + R), P 0,3 m : Percolation below 0.3 m.

Components of Water Balance and Percolation Below 0.3 m Depth
Errors in mass balance (BE) calculated as (Inputs + Outputs + SS) in simulations performed with HYDRUS-2D were minimal, the highest one being 12 mm or less than 1.7% (percentage defined as [ BE Irr+R ·100]). These values are slightly lower than those published by Li et al. [17], where the HYDRUS-1D model was used in rice irrigated by flooding. In general, water balance errors were higher in simulations with silty-clay textured soil.
Cumulative ET c computed following Allen et al. [20] (ET cc ) throughout 240 days of crop cycle accounted 642 mm (Figure 3), while the cumulative ET c simulated with HYDRUS-2D (ET sim ) was among 66 and 86% of the previous one. ET sim was higher in silty-clay and loam soil textures than in sandy-loam and, in general, declined when frequency of between consecutive irrigations decreased and when emitter depth increased, especially in the sandy-loam soil texture. The higher ET sim in loam soil compared with sandy-loam soil texture (cases A and B, Table 4) could be explained by higher grain yield obtained in the part of the field where soil had loam texture. On the other hand, the drained water below 0.3 m depth was higher in the sandy-loam soil than in the loam and silty-clay soil textures, which explains why simulated evapotranspiration and ET sim /(Irr + R) ratios were higher.
In cases where irrigation restored crop evapotranspiration (cases D-V, Table 4) and dripline depth was 0.15 m (D, E, F, J, K, L, P, Q, and R), an increase of the ratio between evapotranspirated water and supplied water ET sim /(Irr + R) was determined in comparison with irrigation patterns observed in the experimental tests (cases A, B, and C), whatever the texture or frequency of irrigation was. This would indicate that following the irrigation criteria to restitute the crop evapotranspiration would have improved the water use by the crop.
Frequency of irrigation had practically no influence on the ET sim in silty-clay soil texture (cases P-V). On the other hand, a slight increase in ET sim /(Irr + R) ratio was observed when increasing the frequency of irrigation in the sandy-loam and loam soil textures. Similarly, an increase in the irrigation frequency reduced drainage below 0.3 m depth. The lower water holding capacity of coarse textured soil may explain the obtained results.
When comparing the 2 dripline depths in a sandy-loam textured soil using the same irrigation frequency, it was observed that 0.15 m depth reduced drainage below 0.30 m, increased the ET sim /(Irr + R) ratio and, consequently, the water used by the plant. This also happened in soils with loamy and silty-clay textures, although to a lesser extent.

Yield and Water Productivity
Yield in the field showed two differentiated zones corresponding to the diverse soil texture within the same field. In the sandy-loam textured area, yield (734 kg ha −1 ) was significantly lower (p < 0.05) than the one corresponding to the loam-textured area (5363 kg ha −1 ). As mentioned in Section 3.1, soil water contents were measured in the sandy-loam part of the field. In the mid-way between two driplines the soil water content achieved values corresponding to matric potential bellow −20 kPa during the period from 34 to 70 DAS (corresponding to from tillering to beginning of flowering stages). Rice is especially sensitive to lack of water during this period, explaining the low yield obtained in the parts of the plot of sandy-loam texture. Furthermore, low water availability in the mid-way between two consecutive driplines in sandy-loam-soil (position 3 in Figure 1), could give a competitive advantage to weeds that have the capacity to develop despite having low amounts of water available. However, the yield of 5363 kg ha −1 and the water productivity (WP Irr+R ) of 0.62 kg m −3 in the loam textured part of the field was similar to that obtained in previous irrigation campaign in a field test irrigated by surface drip in a silty-clay textured soil with a very shallow aquifer (Table 5), in which a yield of 5565 kg ha −1 and a WP Irr+R of 0.60 kg m −3 were obtained [10]. Table 5 shows the water productivity in the current field tests with SDI and the results obtained in continuous flooding irrigation (CFI) and surface drip irrigation during the 2017 campaign in the same area.
Under continuous flooding irrigation (CFI) and dry seeding conditions in nearby fields, the same authors reported a yield of 6486 kg ha −1 and a WP Irr+R of 0.42 kg m −3 . Although in CFI the yield was 16% higher than in SDI, the WP was similar to that of surface drip irrigation (DI), showing an increase above 40% of the WP Irr+R compared with CFI. Therefore, SDI systems may be of interest in situations where water is a limiting factor. Table 5. Grain yield, irrigation water productivity (WP Irr ) and water productivity considering irrigation and rain (WP Irr+R ) in subsurface drip irrigation (SDI), surface drip irrigation (DI) and continuous flooding irrigation (CFI) in Baix Ter area. Among the few published works of rice irrigated using SDI, Rajwade et al. [13] in a sandy-loam soil obtained yields ranging from 1881 to 5074 kg ha −1 and WP Irr+R from 0.23 to 0.61 kg m −3 . Parthasarathi et al. [9] in a clay-loamy soil obtained yields of 5600 kg ha −1 with WP Irr+R of 0.86 kg m −3 . Although the maximum yields reported by these authors were very similar to each other and to the one obtained in the present work, there is a great difference among the WP values, which could be attributed to different soil textures. In general, coarser soils require higher irrigation water amounts in order to allow water to reach the root system from the emitter.

Conclusions
To maximize potential water savings of subsurface drip irrigation (SDI) systems, it is necessary to optimize its design and irrigation management, highlighting depth of emitters and frequency of irrigation. The RMSE and R 2 values obtained by comparing the water contents measured in the field and simulated with HYDRUS-2D showed a reasonably good prediction with the latter, indicating that numerical simulations can be: (1) a useful tool to optimize the design and management of SDI taking into account the hydraulic properties of soils, as well as, (2) predicting water percolation. The results of the simulations applied to the production of rice in Baix Ter area considering soils with different texture allowed to determine that the most adequate dripline depth to minimize water losses due to percolation and maximize water extraction by plants is 0.15 m with a high frequency of irrigation, such as one or two irrigation events per day.
On the other hand, location of measurement points is critical when using irrigation criteria based on a matric potential threshold. Indeed, irrigation management based on maintaining a certain potential or water content in a very close location to the emitter can lead to insufficient irrigation. For this purpose, simulations with HYDRUS-2D can be useful to determine the most appropriate location of the soil water probes in order to efficiently manage the SDI in rice.