Using Geomatic Techniques to Estimate Volume–Area Relationships of Watering Ponds

: Watering ponds represent an important part of the hydrological resources in some water-limited environments. Knowledge about their storage capacity and geometrical characteristics is crucial for a better understanding and management of water resources in the context of climate change. In this study, the suitability of different geomatic approaches to model watering pond geometry and estimate pond-speciﬁc and generalized volume–area–height (V–A–h) relationships was tested. Terrestrial structure-from-motion and multi-view-stereo photogrammetry (SfM-MVS), terrestrial laser scanner (TLS), laser-imaging detection and ranging (LIDAR), and aerial SfM-MVS were tested for the emerged terrain, while the global navigation satellite system (GNSS) was used to survey the submerged terrain and to test the resulting digital elevation models (DEMs). The combined use of terrestrial SfM-MVS and GNSS produced accurate DEMs of the ponds that resulted in an average error of 1.19% in the maximum volume estimation, comparable to that obtained by the TLS+GNSS approach (3.27%). From these DEMs, power and quadratic functions were used to express pond-speciﬁc and generalized V–A–h relationships and checked for accuracy. The results revealed that quadratic functions ﬁt the data particularly well (R 2 ≥ 0.995 and NRMSE < 2.25%) and can therefore be reliably used as simple geometric models of watering ponds in hydrological simulation studies. Finally, a generalized V–A power relationship was obtained. This relationship may be a valuable tool to estimate the storage capacity of other watering ponds in comparable areas in a context of data scarcity. generalized V–A relationship was obtained from the surveys carried out at eight small watering ponds from three representative rangeland farms in the SW Iberian Peninsula that may be used to estimate the storage capacity of other watering ponds in this region.


Introduction
Watering ponds are the main source of drinking water for livestock in the rangelands of the southwestern (SW) Iberian Peninsula [1]. Most of these ponds consist of small earth dams that collect surface runoff from ephemeral streams, with pond sizes rarely exceeding 1 ha [2,3]. Understanding the hydrological functioning of this type of infrastructure is crucial for efficient water management in extensive livestock farms, especially in semi-arid areas, where water resources are often scarce.
Simulation of water-level fluctuation through water-balance models is a frequently used approach to analyze the hydrological dynamics of small water bodies [4]. In these models, the temporal variation of both the stored volume (V) and the flooded area (A) in the water body are often required inputs. However, due to the difficulty of directly monitoring V and A, these are usually calculated from the height of water surface above the pond bottom (h), using predetermined V-h and A-h relationships [5]. It is usual to define these relationships using simple geometric models. The most frequent models use power or quadratic functions to express the V-h and A-h relationships, which has been successfully tested in simulation studies of reservoirs [6], wetlands [7,8], and lakes [9]. V-A relationships (usually defined by power functions) also can be used to model the morphometry of individual water bodies in hydrological simulation studies [10][11][12][13]. In fact, this is the method implemented in the widely used SWAT model [14] to morphologically define reservoir-type elements.
V-A-h relationships can be developed for a particular pond if a detailed bathymetry is available [4]. The most common methods traditionally used to survey the bathymetry of small water bodies include total station, the global navigation satellite system (GNSS), and meter stick [15]. However, the use of these methods may be too costly, labor intensive, and time consuming, especially when a large number of water bodies have to be surveyed.
Laser-based survey techniques evolved quickly in the last decades, and they have been widely used to provide topographic data for hydrological applications [16]. At this point, we may differentiate between aerial and terrestrial systems, known as laser-imaging detection and ranging (LIDAR) and terrestrial laser scanner (TLS), respectively. There are freely available national LIDAR datasets (e.g., the National Geographic Information Center in Spain [17]), but their use to produce reliable V-A-h relationships that can be used in hydrological studies is limited [8,18]. Terrestrial systems (fixed or mobile) may capture terrain characteristics with higher resolution and accuracy than aerial sensors [16,19], but the cost and operation of these instruments are important limiting factors [15,19,20].
On the other hand, the characteristics of livestock watering ponds in the SW Iberian Peninsula (size, water depth, and turbidity) make unsuitable the use of technologies commonly used to survey submerged areas, such as bathymetric green airborne LIDAR [21], bathymetric sonar on manned boat [22], or through-water TLS [23]. Nonetheless, the rapid development in recent years of low-cost mobile platforms such as unmanned aerial vehicles (UAVs) and unmanned surface vehicles (USVs), with increasing use in surveying applications, is making it possible to solve many of the drawbacks associated with traditional surveying techniques in the presence of shallow and turbid waters. For example, a small unmanned vessel can be equipped with a bathymetric sonar to survey the submerged terrain [24], while the emerged areas can be surveyed using UAV photogrammetry [25][26][27][28][29].
Low-cost and efficient methodological alternatives based on satellite data have been proposed by several authors to estimate V-A-h relationships for lakes or reservoirs [30][31][32]. These methods use simultaneous measurements of the flooded area and the water level, the latter being determined in situ [12,33,34] or by means of satellite altimetry [35][36][37]. However, although remote-sensing techniques have proved useful for detecting and mapping small ponds [38], these approaches do not provide sufficient accuracy for topo-bathymetric applications in water bodies with a surface area <1 ha [39][40][41].
Some recent developments in 3D photo-reconstruction techniques, such as the concurrent use of structure-from-motion (SfM; [42]) and multi-view-stereo (MVS; [43]), have contributed to the rapid and cheap production of high-resolution point clouds with similar accuracy to that provided by TLS and conventional photogrammetry [44][45][46][47]. Livestock watering ponds in the SW Iberian Peninsula seem to be suitable for SfM-MVS techniques, as they commonly lack vegetation, and the survey may be carried out at the end of the summer, when they are shallow or completely dry. Additionally, these ponds are usually very small, facilitating the acquisition of images from the ground. Examples in the literature on the use of SfM-MVS techniques to model the geometry of water bodies and estimate the V-A-h relationships are scarce. To our knowledge, the only example is that of Langhammer et al. [18], who applied SfM-MVS to images acquired using a UAV platform to estimate the V-h and A-h relationships for an abandoned montane reservoir. Hence, the suitability of terrestrial SfM-MVS photogrammetry to model water-pond terrain has not been tested yet.
The existence of water in the ponds, sometimes throughout the year, makes it difficult to acquire complete and accurate terrain models, requiring combined approaches for the exposed and submerged terrain [48]. In the wider context of hydrology and fluvial geomorphology, examples of the use of a single technique to register exposed and submerged terrain are scarce and have shown significant limitations [49].
Small water bodies like livestock watering ponds are not usually regarded as part of the hydrological system by water agencies, but the accumulation of large numbers of these water bodies can impact the hydrology, and can be the cause of water conflicts because they take water from the available stream flow [50,51]. Therefore, information about the storage capacity of these water bodies is essential for decision-making processes regarding planning and management of water resources. In this regard, regional V-A relationships are often used to estimate and monitor the volume of water stored in non-surveyed water bodies based on data of flooded areas, which can be obtained from databases [52], topographic maps [53], high-resolution aerial photographs [54], satellite imagery [50], or LIDAR datasets [55].
For any of the above applications, a frequent approach in studies at the basin or regional scale is to fit a common geometric model from the known bathymetry of a set of reference water bodies and apply the function thus obtained to the rest of the water bodies in a given region, assuming that the model parameters remain approximately constant if the geomorphological context is uniform [10][11][12][13]50,52,[55][56][57][58][59][60].
On the above considerations, the objectives of this paper are: 1.
To apply and compare different geomatic approaches and techniques to model the topography of small watering ponds (terrestrial or close-range SfM-MVS, aerial SfM-MVS, GNSS, LIDAR, and TLS). Specifically, the suitability of terrestrial SfM-MVS photogrammetry was tested, as it could be a low-cost, high-accuracy alternative to laser technologies or more time-consuming GNSS surveys. Tips on the use of this approach are also provided; 2.
To assess the overall suitability of power and quadratic functions to describe wateringpond geometry by means of pond-specific V-A-h relationships. These relationships could be a valuable tool to be used as a geometric model of watering ponds in hydrological simulation studies; 3.
To obtain a generalized V-A relationship from the surveys carried out at eight small watering ponds that may be used to estimate the storage capacity of other watering ponds in similar rangeland areas.
The paper has been organized in five sections. The first section presents the motivation and objectives of the study. In Section 2, the methodology is presented, including a description of the study area, the field surveys, and the methods used to derive the V-A-h relationships of the ponds. In Section 3, the results of the geomatic measurements are presented, as well as the resulting V-A-h relationships. Finally, Sections 4 and 5 include the discussion of the results and the main conclusions, respectively.

Study Area
For this work, two privately owned farms (Parapuños and La Brava) and a communal farm (Dehesa Boyal Monroy) located in the SW Iberian Peninsula were selected as study areas ( Figure 1). They are representative of the dehesa land-use system, an agro-silvopastoral system characterized by the presence of open wooded pasturelands of evergreen oaks (Quercus ilex subsp. ballota and Q. suber). These rangelands dominate the landscape of SW Iberia (Figure 1a) and are considered of high natural value in Europe, due to the wide range of ecosystem services they provide [61].
The study farms share the dominant geomorphological characteristics of the rangelands in the SW Iberian Peninsula: the landscape is characterized by gently undulating erosion surfaces, incised by small channels with ephemeral flow, giving rise to increasing slope gradients as approaching the main rivers; with slates and greywackes being the predominant bedrocks. The soils are generally shallow, prevailing Cambisols and Leptosols [62]. The climate is Mediterranean, with a humid season from October to May and a pronounced dry and hot season (June-September), particularly in July and August. Rainfall shows a high temporal variability, both annually and inter-annually [63]. Livestock rearing is the main land use in the study farms.

Surveying Watering-Pond Geometry
Throughout the study farms, several ponds were selected in which the water depth was low enough to allow point collection by wading: P1, P2, P3, and P4 at Parapuños farm; B1, B2, and B3 at La Brava farm; and DBM at Dehesa Boyal Monroy farm (Figure 1c-e).
Several geomatic techniques were applied and combined to model the geometry of the selected watering ponds ( Table 1). As stated before, watering ponds present areas that usually remain submerged even at the end of the hydrological year, requiring the support of GNSS or any other technique to survey the inundated areas. In most cases, the submerged areas show an almost flat and horizontal surface that would be easily modeled by surveying a few points. Specifically, terrestrial or close-range SfM-MVS photogrammetry (supported on GNSS for the submerged terrain) was used for ponds B1 and B2; LIDAR and GNSS data were used for ponds P1, P2, P3, and B3; and aerial SfM-MVS photogrammetry (supported on GNSS for the submerged areas) was used for pond DBM. Finally, pond P4 was surveyed by combining TLS and GNSS for the exposed and submerged terrain, respectively. These techniques were used to survey the topography and produce a DEM of each watering pond. Then, the DEMs were tested using independent check control points (CCPs) previously surveyed with a GNSS RTK device. Well-known and state-of-the-art error metrics were calculated comparing the CCPs and the DEMs: the root mean square error (RMSE), the mean error (ME), and the standard deviation of error (STDE). Field surveys of ponds B1 and B2 were carried out in late summer 2019, when the ponds were at their lowest water level (just a few centimeters; Figure 2a). At this time, a GNSS antenna tied to a survey rod could be used to survey the topography of the submerged surface ( Figure 2). Two Emlid Reach RS antennas were operated as base and rover, respectively. The fixed base station registered coordinates from satellites and sent corrections through a long-range (LoRa) link to the rover antenna, which performed in real-time kinematic (RTK) and fixed-solution status modes. Twenty artificially marked additional points ( Figure 2b) were recorded in the emerged area, also by using the rover GNSS in the RTK mode. These points performed later as ground control points (GCPs, which were used to scale and georeference the model and refine camera calibration parameters, n = 15) and CCPs (i.e., check control points; these points were not used in the composition of the model, but rather to check the geometrical accuracy of the 3D model, n = 5) for the SfM-MVS processing. In order to survey the emerged surface in detail, photographs were acquired around the pond with a Canon 550D digital single-lens reflex camera (DSLR), drawing a convergent image network geometry ( Figure 3). The Canon sensor had 5184 × 3456 pixels, and the focal length was fixed to 35 mm. The Pix4Dmapper pro software (v. 4.5.6) was used for the photogrammetric processing with parameters and specifications shown in Table 2. The GCPs were marked after the initial alignment, then the model was reoptimized to support the refinement of camera parameters. At this point, the layer of water was masked from the point cloud. The resulting point cloud for the emerged terrain (Figure 4a-d) was merged with the GNSS data acquired over the submerged surface, and the resulting point cloud was used to produce a DEM of each pond. This task was carried out with the CloudCompare software package and the rasterize tool. This tool allows an estimation of the Z-value for a pixel (i.e., a XY location) that contains several Z-values (which is typical in high-density regions of a point cloud) and, at the same time, allows the interpolation of values in low point-density regions. A pixel size of 0.2 m was defined for the resulting DEMs.  Watering ponds P1, P2, P3, and B3 were surveyed using the GNSS in RTK mode for the submerged surface, while the emerged topography was obtained from the LIDAR dataset provided by the Spanish National Geographic Information Centre [17] (e.g., Figure 5). This dataset shows a point density of 0.5 pts·m −2 and was acquired using a Leica ALS50 sensor. Both point datasets were merged together and then the topographic surface was modelled using the interpolation algorithm topo to raster [65] within the ArcGIS software (10.5), producing a DEM for each pond with a pixel size of 0.5 m (Figure 4b). Additionally, five points were recorded with the GNSS RTK system and performed as CCPs in order to evaluate the accuracy of the resulting DEMs (i.e., these points were not used in the production of the DEM). This quantification includes the error of the instrument techniques (i.e., LIDAR and GNSS) and the interpolation technique (i.e., the topo to raster algorithm).
Watering pond DBM was surveyed using a fixed-wing UAV (Ebee classic by Sensefly) with a Sony WX220 sensor on board (18 Mpx) in early summer. At this time, a water layer of a couple of decimeters covered the bottom of the pond. The topography of DBM was obtained using a dataset previously acquired by Alfonso-Torreño et al. [66] for the exposed terrain and, again, the submerged terrain was surveyed by means of the GNSS RTK system. A total of 84 points were collected by means of the GNSS, from which 79 were used in combination with the point cloud derived from SfM-MVS to generate a DEM, with the remaining 5 points employed as CCPs. A pixel size of 0.2 m was selected for the resulting DEM.  Finally, watering pond P4 was surveyed combining a TLS (Faro Focus 3D X330) and the GNSS in RTK mode for the emerged and submerged areas, respectively. A total of 15 stations were registered together using a set of artificial spheres (diameter of 14.5 cm). In order to georeference the TLS data, 8 artificial points were marked and surveyed with the GNSS. Finally, 706 points were surveyed in the submerged area by means of the GNSS. Ten additional locations were surveyed in the study area to be used as CCPs to evaluate the quality of the resulting DEM. This DEM was produced by merging together the TLS cloud and the GNSS data, using the rasterize tool within CloudCompare software and defining an output pixel size of 0.2 m.
Maximum uncertainties in the estimation of the water-storage capacity of the ponds were calculated from the maximum area of each pond and the vertical error estimated for each approach, using the following equation: where V E is the error associated to the estimation of the water storage capacity of a pond, RMSE z is the root mean square error of the Z coordinate calculated for the method used to model the topography, and A max is the maximum area of the pond. Specifically, the RMSE z was calculated using the CCPs surveyed by the GNSS RTK system, while A max was estimated using the DEM and the altitude of the spillway (overflow channel).

Obtaining Volume-Area-Height Relationships
For each of the study ponds, V and A values were calculated at intervals of 10 cm in height (over the entire depth range) from the DEM of each pond using the geoprocessing tools of the ArcGIS v10.5 (ESRI) software. In each case, the values of V and A calculated at the maximum water height (i.e., at the elevation of the pond spillway) were assigned the maximum volume (V max ) and the maximum surface area (A max ) of the pond, respectively. Then, Microsoft Excel software was used to derive pond-specific V-h, A-h and V-A relationships by fitting (by least squares) power and quadratic functions of the forms: with α, β, a, b, and c being the adjustment parameters. The goodness of fit of the resulting relationships was determined with the coefficient of determination (R 2 ) and the normalized root mean square error (NRMSE). The latter was calculated for volume (NRMSE V ) or area (NRMSE A ) as follows: where i is the index number for the pond water level, k is the total number of pond water levels, V DEM is the pond volume derived from the DEM at water level i, V F is the pond volume derived from the fitted function (either the power or the quadratic) at water level i, A DEM is the pond area derived from the DEM at water level i, and A F is the pond area derived from the fitted function at water level i. In order to obtain generalized V-h, A-h, and V-A relationships that may be used at the regional scale, the values of V, A, and h (obtained every 10 cm depth) from all study ponds were represented together and fitted equally by means of power and quadratic functions.
Finally, a generalized V-A relationship was also derived by fitting a power function to the V max and A max values of the study ponds, to test whether sufficiently reliable V max values could be obtained based solely on the A max values of the study ponds.

Suitability of Terrestrial and Aerial SfM-MVS Photogrammetry to Model the Topography of Small Watering Ponds
The point clouds produced for the emerged parts of ponds B1 and B2 showed huge amounts of points and resulted in point densities of 19,143 and 24,515 pts·m −3 for B1 and B2, respectively. These volumetric point densities allowed the identification of microtopographic features and very fine details (Figure 4a-d and Table 3). The point density was uniform in the area of interest as the result of a well-planned and executed convergent image network geometry (Figure 3). The error statistics estimated for B1 and B2 showed an outstanding performance of the technique, with RMSEs always below 0.03 m (Table 3). On the other hand, the aerial SfM-MVS produced a lower density (39 pts·m −3 ) and a less-accurate (RMSE CCP = 0.092 m) point cloud (Figure 4g,h). Table 3. Characteristics of 3D point clouds produced by means of SfM-MVS photogrammetry (ponds B1, B2, and DBM) and the TLS instrument (pond P4). GSD: ground-sampling distance; GCP: ground control point; CCP: check control point; RMSE: root mean square error; ME: mean error; STD: standard deviation of error for CCP; UAV: unmanned aerial vehicle. * Note that data for DBM refer to a larger study area and dataset used in Alfonso-Torreño et al. [66].

B1
B2 DBM * P4 The concurrent use of SfM-MVS and GNSS (for the submerged terrain) produced very accurate DEMs for the terrestrial datasets. The estimated RMSE CCP-z for these DEMs resulted in an average error of 1.19% in the maximum volume estimation, varying from 0.72 to 1.65%. The RMSE CCP-z for the DEM elaborated using the aerial SfM-MVS dataset and the GNSS was an order of magnitude larger (RMSE CCP-z = 0.128 m) than the terrestrial approach, and resulted in an error of 12.15% in the estimation of the pond capacity. Note that the RMSE CCP shown in Table 3 are a final independent validation of the resulting DEM surfaces. Therefore, this parameter included the interpolation error and the error associated to any technique used to produce the DEM (i.e., SfM-MVS photogrammetry and GNSS).

GNSS, TLS, and LIDAR
The RMSE of the locations surveyed by the GNSS at the submerged surface of the watering ponds varied from 0.007 m to 0.017 m (Table 4). These were expected figures for a widely used and tested surveying technique, and differences between ponds were not statistically significant. No spatial pattern was observed in the RMSE of the GNSS-surveyed locations, suggesting a uniform spatial distribution of errors. Note that the RMSE figures shown in Table 4 refer exclusively to the accuracy of the GNSS locations registered, while the interpolation errors of each DEM are shown as RMSE z in Table 5. This parameter was calculated using independent CCPs (acquired also by means of the GNSS) and the DEMs that were used later to estimate V max and A max ( Table 5).
The TLS produced by far the densest and most accurate data among the techniques and instruments used ( Table 3). The average point density of the TLS-derived point cloud was 58,366 pts·m −3 ; i.e., more than twice the estimated density of point clouds produced by the terrestrial SfM-MVS techniques. However, point densities were highly variable spatially, with larger values close to the TLS stations (Figure 4e,f). The average registration error of the different point clouds-stations of the TLS instrument was 0.002 m, varying from 0.0003 m to 0.008 m. However, the need to merge (and georeference) the TLS-acquired point cloud with the GNSS data recorded for the submerged terrain produced a degradation of this accuracy to the GNSS range of errors (0.01-0.03 m). The DEM produced for pond P4 with the combination of TLS and GNSS showed a RMSE z of 0.036 m that resulted in an error of 3.27% in the estimation of the maximum volume. The watering ponds modelled using LIDAR and GNSS techniques showed a DEM with an RMSE z that varied from 0.109 m to 0.333 m. These figures represented and average error of 17.51% in the estimation of the V max , varying from 10.86% for pond P2 to 31.08% for pond P3. Table 6 and Figure 6 show the pond-specific V-h, A-h, and V-A relationships that resulted for the study ponds when power and quadratic functions were fitted. In general, good fits could be obtained with both power and quadratic functions (R 2 > 0.9), indicating that watering pond morphometry was adequately described by these types of functions. However, a significantly better performance of the quadratic functions was evidenced, especially for the V-h and A-h relationships. Thus, while a NRMSE greater than 10% was obtained in some ponds when the potential functions were used, this statistic rarely exceeded 2% with the quadratic functions (Table 6). These low errors were accompanied in all cases by excellent coefficients of determination (R 2 > 0.995). Power functions gave rise to fits with such levels of goodness only in the case of V-A relationships, with R 2 values higher than 0.99 and NRMSE around 5% (except for pond B3: R 2 = 0.985, NRMSE = 8.54%). Of the eight watering ponds analyzed, six (75%) had power V-h relationships in which the parameter β was less than 3 ( Table 6), indicating that watering ponds tend to have concave bathymetries [67]. Also noteworthy, in the case of power functions, was the similarity of the adjustment parameters between ponds of the same farm (Table 6).

Pond-Specific and Generalized V-A-h Models
Regarding the generalized V-h, A-h, and V-A relationships (Table 7 and Figure 7), the resulting coefficients of determination and the NRMSE values were of similar magnitude for power and quadratic functions. In general, the fits were quite poor in the case of the generalized V-h and A-h relationships. However, relatively good fits could be obtained for the generalized V-A relationships, especially when the quadratic function was used (R 2 = 0.863 and NRMSE values ranging between 5% and 18%; Table 7).    Figure 8a presents the generalized V max -A max relationship derived from the V max and A max values of the study ponds, expressed by the following equation:

Generalized V-A Relationships for Estimating Water-Storage Capacity
where V max and A max are expressed in m 3 and m 2 , respectively. This function fitted the data reasonably well (R 2 = 0.83), although the volumetric errors, when used to estimate the maximum capacity of the ponds, were highly variable (Table 8), ranging from −2.4% (pond DBM) to 34.9% (pond P1), and being the average of the absolute values of 17.7%. In any case, these estimates were more accurate than those obtained using the generalized V-A relationships derived from the complete geometry of the ponds, whose equations are shown in Table 7. With the latter, whether adjusted by using power functions or quadratic functions, the errors associated with estimating the maximum capacity of the ponds reached values above 40% in many cases (Table 8).  Table 8. Maximum pond capacities (V max ) estimated using the generalized V-A (power and quadratic) relationships, the generalized V max -A max relationship, and the farm-specific V max -A max relationships; and associated volumetric errors (V E ) in relation to the V max values obtained from the DEMs. Interestingly, very good fits were found at the farm level (R 2 > 0.96 and |V E | < 8%; Figure 8b,c and Table 8), which suggested that some factors acting at that scale were relevant in determining the morphometry of the ponds.

Discussion
The combination of terrestrial SfM-MVS photogrammetry (for the emerged terrain) and GNSS (for the submerged terrain) probed to produce accurate topographic models of small watering ponds, with errors in the range of 1-3 cm. These findings were in agreement with other specific previous applications of the SfM-MVS technique to reconstruct similar features (topographic depressions), such as channels or gully headcuts, also using convergent image network geometries [68,69].
The use of the GNSS to support and to complement the SfM-MVS technique was necessary in our study areas due to several reasons. The first was that, even at the end of the summer, watering ponds usually show a thick layer of water at the bottom, which is usually turbid [70] and makes unsuitable underwater photo-reconstructions such as those carried out successfully in other environments [71]. The second reason was that the photogrammetric model needed to be (at least) scaled using GCPs (or GNSS RTK systems connected to the camera instead of GCPs; e.g., [72]), while the GNSS allowed the survey of additional CCPs to independently quantify the accuracy of the model. In recent years, affordable (or low-cost) GNSS survey-grade devices have been available in the market [73] to carry out this task. For example, the two antennas (rover and base station) used here (Emlid Reach RS, Figure 3a) have a price of ≈ EUR 800 (each), and they may be controlled with any smartphone or tablet device. The last reason was that, even without water, the incidence angle that resulted from the hand-held-camera acquisition approach and the almost flat topography of the bottom of watering ponds (Figure 2a) may have been glancing, reducing the number and density of tie points [74] and influencing the accuracy of marking GCPs and CCPs [75]. The use of the GNSS to cover the submerged terrain is well known in the literature and produces accurate bathymetric models [76]. However, GNSS surveys are time-consuming, and point densities are usually lower than those of laser-based or photogrammetric technologies, as evidenced in our study.
The development of SfM-MVS photogrammetry in the last decade has democratized the use and production of high-resolution DEMs [77]. The present application shows that models produced by this technique may provide accurate estimations of pond topography and storage capacity. Once the DEM is produced and the V-A-h relationships are determined for a specific pond, a permanent surveying rod would be sufficient to know the volume of water stored in almost real time. The availability of this data may lead to a better farm management, particularly in water-limited environments such us dehesas [3]. The use of terrestrial SfM-MVS and GNSS, instead of LIDAR and GNSS, represents an important improvement, as it reduces errors in an order of magnitude, allowing the accurate estimation of water stored in the ponds even in periods with low levels.
Our results indicated that the use of existing LIDAR datasets in combination with GNSS only provided rough estimations of water-storage capacity and were particularly unsuitable in periods of water scarcity, when errors were larger than the amount of water stored.
In the complete absence of water in the pond, the use of UAV platforms in combination with SfM-MVS photogrammetry may become an effective alternative that overcomes the glancing perspective of the hand-held camera in the terrestrial approach [74]. In this sense, UAVs that integrate RTK or post-processing kinematic (PPK) systems on board are a promising technology that may reduce the time-consuming and labor-intensive task of GCP-CCP acquisition [78]. Any approach supported on GNSS data acquisition should balance number, spatial extent, and point density. For instance, GNSS devices are suitable to model the flat and horizontal bottom of the ponds, but unsuitable to survey steep and rough banks.
The TLS is the instrument that produces the most accurate and densest point clouds; however, the need to register the submerged terrain or to georeference the resulting models will degrade the accuracy to that of the georeferencing technique or instruments [79]; for example, GNSS systems. Experiences of shallow-water bathymetry in conjunction with emerged ground surveying by means of TLS instruments are scarce and limited to specific hydraulic and water-quality conditions [49], which is not the case for watering ponds in the SW Iberian Peninsula.
On the other hand, the resulting TLS point clouds commonly show heterogeneous point densities (larger close to stations) and occlusions due to the line-of-sight effect [69]. The TLS coverage may be improved by increasing the number of TLS locations [69], which is not always efficient, as the setup of the equipment and data acquisition may take approximately 30 min for each station, and the resulting cloud may show unnecessarily high point densities in other regions. Additionally, the cost of the equipment is still an important limitation for most users [79]. The camera-based SfM-MVS method provides better coverage in very complex topography or rough surfaces [69,80]. On the other hand, there are some experiences that open the door to the use of a single TLS to survey submerged and emerged terrain, but only under specific hydraulic and physical waterquality conditions [49], and this is not the case for watering ponds in the study area.

Suitability of V-A-h Models to Describe Watering Pond Morphometry
This study demonstrated that simple analytical models based on quadratic or power functions can suitably describe watering-pond morphometry. The low associated errors allow these models to be used in hydrological simulation studies. Nonetheless, for simulation models using V-h and A-h relationships, we recommend the fits to be carried out preferably with quadratic functions, for which the NRMSE was less than 2% in almost all watering ponds in our study (Table 6). If V-A relationships are used, both power and quadratic functions would be appropriate (NRMSE rarely exceeded 6% in either case). The high coefficients of determination (R 2 ) obtained for the pond-specific V-A-h models were in line with those reported for small reservoirs in other regions [81,82].
In an attempt to develop a V-A-h model that could be applied to non-surveyed watering ponds in the region, generalized V-h, A-h, and V-A relationships were derived from the bathymetries of the studied ponds. In general, the high NRMSE values obtained for these relationships (up to 42%; Table 7) evidenced the inconvenience of using regional V-A-h relationships as geometric models in hydrological simulation studies of individual watering ponds. However, relatively high R 2 values were obtained for the generalized V-A relationships (R 2 > 0.8; Table 7), which offers the possibility of using such relationships as geometric models of non-surveyed ponds in basin-wide studies, where the individual errors in the regionalization tend to cancel out at the basin scale [83].
In a previous study carried out in the SW Iberian Peninsula, Marín-Comitre et al. [3] obtained a generalized V-A relationship for livestock watering ponds in the region, expressed by the following power function: where V and A are expressed in m 3 and m 2 , respectively. The coefficients of this equation are similar to those obtained in this study for the generalized (power) V-A relationship (Table 7). However, it should be noted that Equation (7) was derived using only three watering ponds from two different farms. Therefore, the new V-A relationships obtained here, based on eight watering ponds from three different farms, represent a significant improvement compared to that expressed by Equation (7).

Performance of V-A Relationships for Estimating Water Storage Capacity in Watering Ponds
It is known that the V-A relationships of small reservoirs are region-specific and that they vary with the geomorphology and geology of the area [13,50,51]. The local nature of these relationships was confirmed by our results, since similar pond-specific V-A relationships were found between ponds of the same farm (Table 6), while significant differences were found between the equations obtained for each of the farm-specific V max -A max relationships (Figure 8b,c). However, these last differences were greater than expected considering the similarity of the geomorphological and geological conditions between the different study areas, which suggested that other factors acting at the farm level were relevant in determining the morphometry of the ponds. For example, the pond construction methodology could play a decisive role.
Nonetheless, we investigated whether a generalized V-A relationship obtained using several reference watering ponds from various locations could be reliably used to estimate the storage capacity of non-surveyed ponds at the regional scale. For this purpose, V-A relationships can be constructed either using only the values of maximum volume (V max ) and maximum flooded area (A max ) from the reference water bodies [10,50,56,67,84,85] or using the volume and area values over the entire depth range of the water bodies [3,12,13,81]. After testing both methods, it was found that the relationship derived from the V max and A max data (Equation (6)) led to lower errors in the estimation of the pond capacities (Table 8) than the generalized V-A relationships (shown in Table 7), even though similar coefficients of determination were obtained in all cases (around 0.8).
Therefore, we propose the V max -A max relationship expressed by Equation (6) to be used to estimate the storage capacity of watering ponds at the regional level. The R 2 value of 0.83 was comparable to those reported by previous studies for generalized V-A relationships in other regions [11,50,56] and gave, to some extent, confidence in its use, especially for watering ponds lacking other available information. The average volumetric error associated with these estimates amounted to 17.7%, with maximum values above 30% for some ponds (Table 8). This means that the pond capacities estimated from the flooded areas might be inaccurate when those are considered individually. However, it is expected that the individual errors tend to cancel out when the relationships are used for broad-based applications carried out at the basin or regional scale [83], in which a large number of ponds are involved. Ultimately, the proposed relationship can be a useful tool for decision-making procedures when managing water resources in rangeland farms of the SW Iberian Peninsula.

Conclusions
The combination of terrestrial SfM-MVS photogrammetry and GNSS produced accurate topographic models of small watering ponds that were useful to improve existing V-A relationships. These models were suitable to accurately estimate V max , comparable to those produced using TLS+GNSS. The errors of the DEMs and V max estimations produced by the rest of the techniques were, at least, an order of magnitude larger.
Pond-specific V-A-h relationships based on quadratic or power functions can adequately describe watering pond morphometry and may be used reliably as geometric models of ponds in hydrological simulation studies.
A generalized V-A relationship was obtained from the surveys carried out at eight small watering ponds from three representative rangeland farms in the SW Iberian Peninsula that may be used to estimate the storage capacity of other watering ponds in this region.
Funding: This publication has been made possible thanks to funding granted by the Consejería de Economía, Ciencia y Agenda Digital from Junta de Extremadura, and the European Regional Development Fund of the European Union through the reference grants IB20036, GR18037, and GR18053.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available LIDAR datasets were analyzed in this study. This data can be found here: https://centrodedescargas.cnig.es/CentroDescargas/index.jsp (accessed on 24 July 2021).