The Transboundary Nature of the Allende–Piedras Negras Aquifer Using a Numerical Model Approach

The Allende–Piedras Negras (APN) aquifer is located between the states of Texas (United States [U.S.]) and Coahuila (Mexico). The Rio Grande crosses the aquifer, acting as a natural and political divide between the countries. However, it remains unclear whether the APN aquifer can be considered a truly transboundary aquifer flow system, which would potentially require joint management by two different administrative jurisdictions. The main purpose of this study was to evaluate the transboundary nature of this aquifer. This was achieved by developing a detailed hydrogeological model to analyze the direction of volumetric fluxes within the APN aquifer using Visual MODFLOW. The model simulated a spatially averaged cumulative drawdown of 0.76 m for the entire aquifer over an 18‐year modeling period (2000–2017). The flow convergence zone, previously located below the Rio Grande, has shifted to the U.S. side in most locations, driven by higher pumping rates of the wells located near the river. This shift of the convergence zone from one country to the other means that groundwater recharge from one side flows underneath the river to the other side. This qualifies the APN aquifer as a “transboundary groundwater flow system.” The procedure followed in this study may be applied to other aquifers that straddle the U.S.–Mexico border and may motivate future modeling studies on other poorly studied transboundary aquifers around the world and thereby enable bi‐national aquifer management.


INTRODUCTION
A transboundary aquifer is an aquifer in which groundwater flows across geopolitical boundaries separating two or more different administrative units (Eckstein 2017). If there is a river separating two different countries (such as the Rio Grande separating the United States [U.S.] and Mexico), then the aquifer will be a transboundary aquifer if groundwater can move underneath the river from one country to the other country (Rivera 2015). In other words, according to Rivera (2015), the aquifer has a transboundary nature if the converging groundwater flows (water divide) do not match the political boundary, which would theoretically be located in the middle of the river. Flows within an aquifer are not immutable, however, and therefore the transboundary classification of an aquifer may change over time and scale as groundwater flows are influenced by changes in human groundwater use and/or physical conditions such as drought or climate change (Sanchez et al. 2018a). As long as the aquifer has the capability to move groundwater across the river, which defines a political boundary, the aquifer can potentially be considered as a transboundary groundwater flow system (Rivera 2015).
Recently, a total of 36 potential transboundary aquifers have been identified along the U.S.-Mexico border (Sanchez et al. 2016). Sixteen of these aquifers were identified and characterized as transboundary with a reasonable level of confidence, but only 11 aquifers have been recognized officially as transboundary by both countries. The Allende-Piedras Negras (APN) aquifer between Texas (U.S.) and the state of Coahuila (Mexico) has been identified as transboundary with a reasonable level of confidence (Sanchez et al. 2016;Sanchez et al. 2018b). However, it has not been officially recognized as transboundary by Mexico or the U.S. or at the international level (Internationally Shared Aquifer Resources Management -ISARM). There is uncertainty related to whether the APN aquifer should be classified as a transboundary aquifer because of the limited available data and research on the aquifer on both sides of the border. Previous studies have only focused on the central portion of the aquifer, located on the Mexican side, and have excluded the northern portion, located on the Texas side ( Figure 1) (Castillo Aguiñaga 2000;Boghici 2002;Modelo 2003;Aguilar 2013;CONAGUA 2014).
The APN aquifer supplies 85% of the water needs of the region (see Figure 1). The expected socioeconomic development and the growing reliance on groundwater to fulfill agricultural, cattle, and industrial water demands from beer production and maquiladoras, which are Mexican factories along the U.S.-Mexico border, will put additional pressure on the aquifer, which has already shown signs of unsustainable use on the Mexico side (DOF 2011;CONA-GUA 2014). This presents a challenge for future groundwater management at a regional and eventually a binational level.
According to official studies (Boghici 2002;Lesser-Illades et al. 2008), the APN aquifer is an unconfined system that is highly dependent on surface-groundwater interactions to maintain water levels. Depletion of the APN aquifer water table as regional pumping rates continue to increase (DOF 2011; CON-AGUA 2014) could decrease the volumetric flux of groundwater discharging to streams of the Rio Grande system and may substantially decrease streamflow. Although most streams in this dry region flow intermittently, several contribute a significant amount of water to the water budget on the Mexico side. Consequently, these streamflow depletions impact the water allocation obligations set in a treaty between the U.S. and Mexico (IBWC 1944), adding an additional legal and political element to the strategic role of the aquifer at the binational level.
The conditions before the year 2000, when much less groundwater pumping occurred than today, are referred to throughout this paper as the pre-2000 conditions. Under the pre-2000 conditions, the flow pathways within the APN aquifer converged along the Rio Grande because they discharged to the river. As groundwater pumping increased in the areas surrounding the river, the flow convergence zone may have shifted toward well fields with higher relative extraction rates (Boghici 2002). According to Rivera (2015), this type of unconfined aquifer, which contains a river acting as a political border, may have little transboundary flow unless the extraction of groundwater is high enough to reduce baseflow to the river. In this case, the groundwater flow system becomes transboundary as its connection with surface water comes to have a regional hydrological impact in the basin.
The main objective of this work was to demonstrate the hydrogeological linkages of the APN aquifer at the transboundary level to support its identification, recognition, and future management efforts at the binational level. This was primarily achieved through the development of a hydrogeological model using the software Visual MODFLOW (VMF). This model was developed to estimate: (1) water fluxes moving across the political border between Coahuila and Texas; (2) water fluxes moving between the Rio Grande and its tributaries and the aquifer, as well as water discharged to the river from the aquifer, and later infiltration back into the aquifer, known as flow-through (Hoehn 1998); and (3) the effects of pumping wells and extreme weather events on variations in groundwater flow patterns in the region.

Site Description
The APN aquifer is located in the southeast portion of the state of Texas in the U.S. and the northern portion of the state of Coahuila in Mexico. The total area covered by this aquifer is 7,024 km 2 , with 5,427 km 2 lying in Mexico and the remaining 1,597 km 2 in the U.S. (Figure 1a-1c). The aquifer boundaries were delineated to follow the distribution of late Neogene and Quaternary deposits on the lower flatlands of the region (Boghici 2002). These are surrounded by older Paleogene and Cretaceous resistant rocks placed as mountain chains and small hills at the northwestern edge of the APN aquifer (Aguilar 2013).
The study area is located within two physiographic provinces: the Eastern Sierra Madre in the west and the Great Plains of North America in the center and east. The aquifer is bordered by mountain chains known as Serrania del Burro and Lomerio Peyotes in Mexico and the Anacacho Mountains in the U.S. (Figure 1b). These mountain chains consist mostly of outcropping Cretaceous rocks that are separated by flat, elongated valleys, reflecting the calcareous nature of the area. The Cretaceous rocks located in the lowlands were covered by alluvial sediments during the late Neogene and Quaternary, creating the present-day plains in the region (Castillo Aguiñaga 2000).
The main surface drainage in the area is provided by the San Antonio, San Rodrigo, and Escondido Rivers and Castaño Creek, which flow from the Serrania del Burro in Mexico into the Rio Grande. Los Morales Creek flows through Maverick County in Texas and into the Rio Grande. Water flow in these drainage systems is intermittent and forms disconnected, stagnant ponds along their courses. Most of the rivers (except the Escondido River, Castaño Creek, and the Rio Grande) dry up completely during the summer. The Escondido River and Castaño Creek have perennial flows, with average discharges of 4 and 2 m 3 /s, respectively (Aguilar 2013) ( Figure 1b).
Most of the precipitation in the region occurs as sporadic thunderstorms. The predominant climate in the study area is semidry to semiarid (Boghici 2002). According to precipitation records for the period 1960-2007, an average of 446 mm of rain falls each year, with most of the rain falling during the wet season from May through September (Aguilar 2013). In the surrounding cities of Allende, Guerrero, Morelos, Nava, Piedras Negras, Villa Union, and Zaragoza (Mexico), the average annual temperature ranges from 20°C to 22°C, whereas in the highlands of Serrania del Burro, it is 18°C (Aguilar 2013). Maximum average monthly temperatures of 30°C are recorded during July and August in the region between Piedras Negras and Ciudad Acuna (Aguilar 2013).
There are no ground-based measurements of evapotranspiration (ET) available for the study area. The JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION potential ET ranges from an annual average of 1,746 mm in the city of Allende to 1,816 mm in the city of Piedras Negras (CONAGUA 2014). ET was estimated for the region between Allende and Piedras Negras as 433 mm/year by CONAGUA (2014) using the Coutagne empirical equation (Coutagne 1954).

Geology and Hydrogeology of the Study Area
The following section describes the geological features of the formations identified and correlated between Mexico and the U.S. in the study area.
The mountains surrounding the APN aquifer consist mainly of carbonate facies deposited in a reef barrier marine environment during the Lower Cretaceous. Primary porosity in these carbonates is almost absent, but secondary porosity is present in the forms of fractures and dissolution cavities. These impermeable rocks are covered by fine sediment layers of the Upper Cretaceous age that locally form the confined underlying aquifers. The flatlands surrounding the APN aquifer are the product of discordant Paleogene and Neogene continental, fluvial deposits and transitional, swamp, and beach deposits. Quaternary fluvial deposits overlie the Cretaceous and Paleogene rocks unconformably ( Figure 2 and Table 1).
The formations that make up the APN aquifer include the Reynosa Formation (Miocene-Pliocene), the Uvalde Gravel (Pliocene-Pleistocene), and the Quaternary deposits (Pleistocene-Holocene). The APN aquifer is identified in Mexico as Quaternary alluvial deposits and Quaternary conglomerates (Figure 2). It is necessary to clarify that Boghici (2002) classifies the Mexico portion of the aquifer as the Reynosa Formation. However, CONAGUA (2014) and the Servicio Geol ogico Mexicano (2008b) identify the formation as Quaternary Alluvial deposits. In the present study, the nomenclature used for the APN aquifer is from Servicio Geol ogico Mexicano (SGM).
A detailed description of the geological formations that make up the APN is attached in the Supporting Information for this study (Supporting Information-01). In the Piedras Negras region of Mexico, the APN aquifer has been recognized as the main exploitable aquifer by Boghici (2002) and CONAGUA (2014). On the U.S. side, this aquifer is made up of Uvalde Gravel and Quaternary deposits. CONAGUA (2014) describes this aquifer as highly permeable, with transmissivity values ranging from 0.0005 m 2 /s (43.2 m 2 /day) to 0.005 m 2 /s (432 m 2 /day).

Socioeconomic Considerations and Importance of the Aquifer
The region overlying the APN aquifer encompasses important border cities and urban areas. On the Mexico side of the border, the largest urban centers in the region are Zaragoza, Morelos, Allende, Villa Union, Nava, El Moral, Guerrero, and Piedras Negras. On the U.S. side, the cities of Quemado, Spofford, Brackettville, and Eagle Pass are the most important urban centers. In 2010, a total of 337,309 people inhabited the border region, according to census records (CENSUS 2010;INEGI 2010). Piedras Negras and Eagle Pass are the largest urban centers in the area, with a combined population of 272,410 (INEGI 2010).
Agriculture and cattle production are the main economic activities on both sides of the border (Rio Grande Regional Water Group 2015). However, coal mining and other industries, such as beer production, are also important on the Mexico side (INEGI 2010). Several factories, including a firearms factory, are located on the U.S. side; one of these companies, O.F. Mossberg & Sons, Inc., is responsible for more than 90% of Mossberg firearms production, in their facility located in Eagle Pass (Miniter 2014). The development of maquiladoras on the Mexico side since the end of the 1990s increased the population by 3.7% per year in the border area to provide a local workforce (Terry 2017). This has placed enormous additional pressure on the area's natural resources, such as sand and gravel within unconsolidated deposits that are extracted from the San Rodrigo River and used as construction materials (Olivera et al. 2018). Building material extraction activities along the San  Table 1 represent the geological units shown in Figure 2.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
Rodrigo River are permitted by CONAGUA, but illegal extraction also occurs. Riparian zones in particular have suffered deforestation and water quality degradation in the area (Olivera et al. 2018). The APN aquifer supplies 85% of the total water needs of the Mexican portion of the region. Of this total, 69.2% is used for agriculture, 17.3% is used for industrial purposes, 4.9% goes to public water supplies, and the remaining 8.9% is used by rural households (CONAGUA 2014). Groundwater pumping for irrigation and mining activities between the cities of Allende and Piedras Negras has lowered the water table, which has the potential to decrease baseflow to the Escondido River (CONAGUA 2011). In Texas, 46% of the water extracted from the APN is used for irrigation, 26% is used for municipal purposes, and 20% is used to raise livestock (Boghici 2002).

METHODS
To quantify water fluxes across the Rio Grande political border between Mexico and Texas and the impacts of water pumping and extreme weather conditions on changes in groundwater flow patterns in the region, it was necessary to construct a groundwater flow model of the region. This model provided estimates of flowpaths beneath the river (Williams 1993;Sophocleous 2002), the evolution of the water table over an 18-year simulation period, and the volumetric fluxes of water that flows from different regions into the Rio Grande and its tributaries.

Definition of Boundaries
The boundaries of the APN aquifer were defined according to geologic maps at 1:25,0000 scale from the Instituto Nacional de Estadistica y Geografia (Servicio Geol ogico Mexicano 2008a, b), as well as Texas geological map (U.S. Geological Survey 2007). The map used to correlate the geological units from both sides of the border was generated by Page et al. (2009). This map was selected because it contains the lithological formations in Texas and Coahuila with Texas nomenclature, making it possible to correlate them across the border. It was also necessary to use the lithological descriptions for geological units, as well as correlations with the equivalent formations on both sides of the border, which were recorded in the geological maps from SGM ( Figure 2). For the APN aquifer model input, the physical boundaries rely on the lithologic differences between the porous Neogene-Quaternary deposits (Uvalde Gravel, Reynosa Formation, and Quaternary deposits) and the surrounding impermeable Cretaceous and Paleogene Formations ( Figure 2).

Isopachs
After the surface geologic units were delineated, drilling reports obtained for the U.S. side were reviewed. These were downloaded from the public database on the Texas Water Development Board webpage (TWDB 2017). Using the lithological descriptions for the Uvalde Gravel and Quaternary deposits in the area from drilling reports, the aquifer thickness was estimated by matching the lithological descriptions in the drilling reports with geological descriptions from the literature.
No publicly available water well drilling reports were available for the Mexico side. Lithological information from four oil wells, however, was obtained from private companies that have recorded the thicknesses of the quaternary deposits. A geological cross section was used to interpolate the thicknesses of the quaternary deposits near the towns of Nava, Zaragoza, and Piedras Negras (Coahuila). As lithological information was lacking for some areas, the thickness of the APN aquifer was inferred from the distribution of water well depths. Deeper wells tended to be located on the periphery of the APN aquifer, where its expected thickness is <1 m. These wells are assumed to pump water from deeper formations and are not screened within the quaternary deposits or the Reynosa Formation within the APN aquifer. Data from the wells were used to categorize the wells based on the geological formation given in the database. The wells with alluvial aquifer identification codes were used to generate the isopachs map of the Quaternary and Neogene deposits.
The aquifer thickness ranges from 1 to 40 m, according to the generated isopachs. The greatest thickness predominates in the center of the basin, in the region of Coahuila, where most of the pumping wells are located. In the Rio Grande area, the thickness of the aquifer ranges from 10 to 25 m. The river does not intersect the aquifer through erosion; therefore, an important assumption is that Quaternary and Neogene deposits below the river bottom limit river interactions with these shallow deposits. For the southern portion of the study area, there was not enough information available to interpolate the aquifer thickness. However, according to information available from the few wells between Castano and Camaron Creeks, the thickness ranges from approximately 1 m to approximately 30 m (see Supporting Information-02).

Evapotranspiration
Estimates of the monthly average ET over the study area were obtained from the NASA Global Land Data Assimilation System (GLDAS) for the period from January 2000 to December 2017. It is not possible to determine ET directly from remote sensing techniques; GLDAS calculates ET on the basis of a water-energy balance (Rodell et al. 2004). This has a temporal resolution of one month and a spatial resolution of 0.25°9 0.25°(27.5 9 27.5 km per cell). Thirteen of these grid cells fall either partially or fully within the area of the APN aquifer. A comparison of the time series for each cell with its adjacent cells showed that the spatial variability is much lower than the interannual variability. Stated another way, the ET within all grid cells responded together to large seasonal changes such that the differences between cells during any given month were small. Therefore, a single monthly ET value for the entire surface of the APN aquifer was calculated as a weighted average of the values for the cells lying within the APN aquifer, and the same value of ET was assumed for the entire study area to simplify the numerical model and minimize computational time. All hydraulic heads and fluxes were modeled on a monthly basis. Extraction of ET values from NASA GLDAS images on an averaged monthly basis yielded an accumulated averaged estimate of 512 mm/year for the entire study area over the 18-year calibration period. The maximum annual ET of 744.9 mm/year occurred in 2004. The minimum annual ET of 309.8 mm/year occurred in 2011 (see Supporting Information-03).
The extinction depth is defined as the depth at which ET from the water table ceases (McDonald and Harbaugh 1988) and is a necessary value in the calculation of water loss to ET in the model. The extinction depth used in the VMF ET input was 1.5 m, as recommended by Shah et al. (2007) for sandy loam with bare soil, agricultural, and grassland covers.

Precipitation and Recharge
Precipitation data were obtained from the Tropical Rainfall Measurement Mission (TRMM) (3B43 Version 7). The TRMM rainfall estimates are calculated from infrared data (Huffman and Bolvin 2013). Monthly average precipitation data were downloaded from the NASA Giovanni visualization tool (Acker and Leptoukh 2007) for the period January 2000-December 2017, at a spatial resolution of 0.25°9 0.25°(27.5 9 27.5 km). As with the satellitebased ET estimates, the spatial variability of precipitation was much lower than the seasonal or interannual variability. Therefore, a single monthly rainfall amount for the entire surface of the APN aquifer was calculated as a weighted average of the values for the cells lying within the APN aquifer, and the same value of precipitation was assumed for the entire study area to simplify the numerical model and minimize computation times.
Visual MODFLOW requires recharge data rather than precipitation data as input to the model. Therefore, it was necessary to estimate the natural recharge in the area from the monthly precipitation values previously obtained from remote sensing. Recharge was calculated using the U.S. Geologic Survey's Soil Conservation Service (SCS) curve number (CN) method (SCS-CN) (NRCS 2004) and the average monthly precipitation across the study area. According to Stewart et al. (2011), the recommended CN range for semiarid regions is 68.2-92.6. An intermediate value of 80 was used to calculate the potential maximum retention (Equation 1) and initial abstraction (Equation 2) to obtain the amount of water that infiltrates during a rainfall event. The potential maximum retention after runoff starts (S) (mm) is defined as follows: where CN is the curve number. The potential maximum retention represents the rain that will not be converted to runoff. The initial abstraction (I a ) is defined as follows: Initial abstraction (mm) is the sum of the interception, infiltration during early parts of the storm, and surface depression storage.
If the average monthly rainfall amount was greater than the initial abstraction, it was assumed that the remaining water ran off and that only an amount of water corresponding to I a infiltrated. If the average monthly rainfall amount was <I a , it was JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION assumed that all of the water infiltrated into the ground as natural recharge with no runoff.
The infiltrated water may replenish the soil moisture, evaporate, transpire, or recharge the aquifer. Indeed, determination of recharge in a semiarid region can be an extremely challenging task, as application of the conventionally used recharge coefficient concept (the percentage of annual recharge over the annual precipitation) usually fails, as noted by Cheng et al. (2017).

River Stage
River stages were obtained from the IBWC (2018) using four gages: two river gages located on the Rio Grande, one on the Escondido River and one on the San Rodrigo River (see Supporting Information-05). For each gaging station, the river stages at 15-min increments and daily total volumetric discharge were available. Only 15-min stage data for the month prior to data retrieval were available. In contrast, daily volumetric discharges for the period 2000-2016 were available. The datasets with available stage and discharge information for every 15 min were used to calibrate the roughness coefficient (n) in Manning's equation (Equation 3). Approximate channel dimensions (channel width, W; river bottom width, b; and bank zones, C 1 and C 2, ) and slope (Slope) were obtained from Google Earth Pro elevation profiles (Google Earth Pro 2018), considering channels with no geometric variations over the simulation period. After calibrating the roughness coefficient (n) for each of the four river gages, we proceeded to calculate the stage for the datasets for which only daily average discharges were available, using Equation (3). As an extra step, the set of discharge values for 2002 (which was most similar in precipitation to 2017) were used to estimate the missing river gage values for 2017.
The VMF river package accepts river stage (a) as a specified head boundary but does not accept river discharge. Therefore, Manning's equation in SI units (Equation 3) was used to calculate the stage (a) along the channels, using interpolated discharges between stations: where Q is the discharge, A is the flow area, R is the wetted perimeter, and Slope is the channel slope. A trapezoidal channel shape was chosen for the application of Manning's equation because this approximated in most places the best channel shape for streams with gentle-gradient riffle (Rosgen 1994).
The variables A and R may be expressed in terms of channel dimensions to calculate ɑ as shown in Equation (4): where the channel bottom width (b) was assumed to be constant, and the bank zones (C 1 , C 2 , and the distance from the bottom of the river to the water level on the river edge), as well as the channel surface width (W), changed depending on ɑ.
The value of n was calibrated using stage and discharge measurements every 15 min. Then, the values of n, Slope, b, C 1 , C 2 , and W were substituted into Equation (4) to calculate the daily stage (a) using the daily average discharges (Q).
To select a monthly representative value for each river gage, it was necessary to perform a statistical analysis of the daily river discharges and calculated stages for the period 2000-2017. Because the data distribution was not normal and the skewness was high for all the stations, the median monthly stage and discharge values were used.

Initial Conditions and Water Level Observations
A period of 18 years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) was selected because of the lack of data from before the year 2000. Thanks to the availability of datasets for several time periods after 2000, it is possible to perform a transient simulation. Data provided by the U.S. Environmental Protection Agency (USEPA) for 86 observation water levels located in Mexico were used to develop datasets as part of a study by Boghici (2002). On the U.S. side, there was only one observation well, near Spofford, Texas, with water level observation data for the year 2000 (TWDB 2017). The water table elevation was contoured from the total of 87 water levels. River stages for the Rio Grande were used as initial heads on the political boundary, for which reports of well water levels are not available for the period 1999-2000. The water table contour map was assumed to represent the initial conditions on January 1, 2000, because of the lack of data for previous years. These initial conditions, representing a period when pumping activities were not extensive, are referred to as "pre-2000" in this paper, and it remains unclear that these initial conditions existed in equilibrium.
Water level values for the years 2006, 2008, 2011, and 2014, available from 222 observation wells screened within the APN aquifer, were used to calibrate the model. Data for 221 of the 222 wells were obtained from a database provided by private companies with information for Mexico, and data for one observation well in the U.S. were downloaded from the TWDB webpage (TWDB 2017). Of the 221 observation wells on the Mexico side, 170 had only one or two water level measurements. The remaining 51 wells had three or more water level measurements, mostly in the years of 2008 and 2014.

Pumping Wells
Pumping well information was obtained from REPDA (Registro P ublico del Agua), Lesser and Associates, and private Mexican companies for the wells located in Mexico. Pumping well information for the U.S. side was retrieved from drilling reports on the TWDB website. Of a total of 799 pumping wells in the study area, 704 are located in Mexico, and 95 are located in the U.S. Several of these pumping wells in Mexico and the U.S. were reported to be active but had no available pumping rate information. For these wells, a pumping rate of 50 m 3 /day was assumed, because this was the predominant known value for wells in the REPDA databases.
Nine of the wells located in Mexico had reported pumping rates that were anomalously high (>5,000 m 3 /day each). This represents only 1% of the pumping wells. These high pumping rates are questionable for several reasons. First, such high pumping rates produce very rapid drops of the cones of depression in the model, drying the cells where the wells are located, as well as the surrounding cells. Second, data suggest that commercial pumps sold in Mexico during the pump operation have an efficient pumping capacity of approximately 3,000 m 3 /day for wells up to 50 m in depth. Therefore, the pumping rates for these nine wells were adjusted to this value. Third, there may be wells on the Mexico side with unreported volumes. Therefore, to simulate a highstress scenario, a pumping rate of 3,000 m 3 /day instead of the median or the mode (50 m 3 /day) was preferred.
The pumping rates for 92% of the wells were close to 50 m 3 /day. The remaining 8% of the wells had rates of more than 750 m 3 /day. In the absence of well screen intervals in the databases, it was assumed that the wells fully penetrate the APN aquifer. The pumping schedules were assigned in the model based on comments provided by institutions managing the well information. If no pumping schedule information was available, it was assumed that a well was pumped constantly at the assigned rate in the REPDA or TWDB database. These assumptions were necessary to complete the information regarding pumping rates and schedules; however, they represent a source of uncertainty in the model.

Layer Definition
According to Castillo Aguiñaga (2000), it is possible to identify two continuous layers in the alluvial deposits between Lomerio Peyotes and Serrania del Burro (Figure 1). In addition, to observe the water flow below the river bottom, evaluate water fluxes moving across the border and in and out of the Rio Grande, and examine the influence of pumping wells, it was necessary to set the river cells only in the upper layer, whereas most of the wells were assumed to be screened in both layers. This distribution in two layers allowed the representation of the rivers and interactions with the shallow zones in the APN aquifer in the upper layer, assigning the lower layer exclusively to groundwater interactions of the APN aquifer below the Rio Grande system (Figure 3).

Hydraulic Conductivity
Boghici (2002) estimated that the horizontal hydraulic conductivity (K x ) ranges from 160 to 430 m/ day near the Quemado Valley adjacent to the Rio Grande. The K x value used in the riparian region of the Rio Grande was 160 m/day, and this same value was also assigned to the riparian areas of the Escondido and San Rodrigo Rivers (Boghici 2002). Castillo Aguiñaga (2000) divided the aquifer into two layers with different K x values that were obtained from pumping tests performed in 1996. A range of K x between 64 and 128 m/day was calculated for the upper layer. In contrast, K x ranged from 7 to 91 m/ day for the lower layer.
In this study, the lowest values of 64 m/day for the remaining upper layer and 7 m/day for the lower layer were selected, as recommended by Hill and Tiedeman (2008), who found that the estimated K x,y values obtained from numerical models are often smaller than the measured values obtained from pumping tests because the gravel pack surrounding the well screen has a higher K than the surrounding media ( Figure 3).
In the APN aquifer, a combination of values from previously cited studies were used to determine distinct K zones. The K value for VMF must be defined in the X, Y, and Z directions. The relationship between K x and K z is known as the hydraulic anisotropy ratio, and a value of 10 is recommended for JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION alluvial deposits with clay lenses (Todd 1980). In the model, it was assumed that K x = K y = 10K z . Vertical anisotropy in alluvial aquifers, however, can range from values of more than 335 on the Ganges Brahmaputra Delta (Knappett et al. 2016), to values closer to 1, measured for the Mississippi River Alluvial aquifer (Carlson 2007).

Specific Yield, Total Porosity, and Effective Porosity
Boghici (2002) recommended a specific yield (S y ) of 0.22. This is a typical value for unconsolidated sandy materials such as the APN aquifer (Heath 1983) and was selected for the numerical model. A porosity (/) value of 0.25 was used, as recommended by Heath (1983) for unconsolidated sands.
The assumption of both total and effective porosities being equal is a common practice in numerical modeling of sandy aquifers (Zheng and Bennett 2002).

Water Budget Zone
Three water budget zones were defined over the study area. One zone was defined for the Texas (U.S.) side, one was defined for the Coahuila (Mexico) side, and one was defined for the Rio Grande corridor. The aim was to evaluate the water exchanges between the river corridor and the regions of the aquifer to the north and south. A time step of 0.1 days near the start of the transient model simulation was used to evaluate water fluxes in the aquifer without pumping. The pumping wells were programmed to start after 0.1 days to ensure they would not have any influence on the pre-2000 conditions in the model, including the water budget before intensive pumping activities began.

Description of the Conceptual Model
Visual MODFLOW was used to generate a numerical groundwater flow model to evaluate flow patterns in the APN aquifer, while primarily focusing on the Rio Grande river system and the surrounding corridor. The study area was 218 km long and 76 km wide, covering the entire APN aquifer. The grid size was set to 500 9 500 m, with 436 rows and 152 columns. The border of the aquifer was set as unbounded. It was assumed that no cross-formational flow occurred from underlying aquifers, because of a lack of information regarding this matter. The aquifer thickness varies from 1 to 40 m and was divided into two layers, as described by Castillo Aguiñaga (2000) for the Reynosa Formation, which is the predominant geological unit in the aquifer. The thickness of the upper layer ranges from 1 m in the western part to 25 m in the Rio Bravo/Rio Grande and includes the upper Reynosa Formation and the Quaternary deposits. The lower layer ranges in thickness from 1 to 15 m and is exclusively composed of Reynosa Formation material. The purpose of this division was to properly represent the river location and the hydraulic parameters described in the Methods section. The main sources of recharge to the APN aquifer are precipitation and infiltration from periodic and perennial rivers. There are areas, however, where the water level increases even during dry periods (low recharge) and far from the influence of perennial rivers and streams. These water level changes could be an indicator of cross-formational flow from the Cretaceous formations underlying the APN aquifer, as suggested by Boghici (2002) and Modelo (2011). This particular observation is discussed in the Results section.
The discharge points from the APN aquifer include pumping wells used in irrigation and municipal supply, ET, and most of the rivers. The Rio Grande is mostly a gaining river in the APN aquifer, which means that the aquifer discharges water into the river. For example, the river gains water from the aquifer near Quemado, Texas and El Moral, Coahuila. In the area north of Piedras Negras, Coahuila, and Eagle Pass, Texas, however, the direction of flow is reversed. Here, the Rio Grande loses water, recharging the aquifer. This reversal is driven by the high density of pumping wells between Piedras Negras and Villa Union (Coahuila), which drive water table declines in the area. An additional loss of water from the APN aquifer occurs through ET, which is significant in this arid to semiarid climate. This eventually reduces the effective recharge. Figure 4 shows the conceptual model of the APN aquifer, including the Texas and Coahuila sides. There is evidence that there is potential cross-formational flow from underlying carbonate aquifers, such as the Austin Chalk and the Edwards-Trinity Aquifer. Batzner (1976) suggested that cross-formational flow from the Edwards-Trinity aquifer is recharging the APN aquifer near the Lomerios Peyotes area. This was also suggested by Castillo Aguiñaga (2000), Boghici (2002), and Lesser-Illades et al. (2008), based on the prevalence of carbonatetype groundwater between Guerrero and Villa Union (Coahuila), in the surrounding areas of Lomerio Peyotes, and in the northern central portion of the APN aquifer ( Figure 4).
The locations with high carbonate content were identified as Type 1 and marked on the figure "Hydrochemical facies distribution in the APN aquifer" in the Supporting Information-06 (Boghici 2002). A similar pattern in the APN aquifer portion located in Mexico was observed by Castillo Aguiñaga (2000). These are labeled as Zones 1 and 2 in a map in the Supporting Information-09 and -10.
It is important to note that the model was run without assuming cross-formational flow, because of lack of information on cross-formational water volumes. Therefore, groundwater volumes recharging the APN aquifer from other aquifers should be quantified for more accurate calculation of modeled fluxes on all boundaries of the APN aquifer.

Numerical Model
To better understand the APN aquifer as a groundwater system, a numerical model was created in VMF to represent the development of the aquifer after the year 2000. A transient model was chosen and was run for 18 years (6,575 days). The selection of the solver to run the numerical model was based on comparisons of the accumulated mass balance at the end of the simulation, as suggested by Kumar (2015). The Waterloo Hydrogeologic Solver indicated a lower cumulative discrepancy (3.7%) than the Preconditioned Conjugate Gradient solver (4.9%) and was therefore selected to run the numerical model.

Calibration
Before calibration of the model, a sensitivity analysis was conducted to test the ranges of K x , /, S y , and recharge. One of the findings of the sensitivity analysis was that varying K x from 0.01 to 1,000 m/day did not change the water table height significantly unless the values were modified by several orders of magnitude. Even then, the water table height only changed by several millimeters, and this may be explained as follows. Hydraulic diffusivity for an unconfined aquifer is defined as a ratio of transmissivity (a product of K and saturated thickness) over the specific yield. A higher value of hydraulic diffusivity implies that the flow transiency will last a short period, and vice versa (Barlow and Leake 2012). The flow system will become less sensitive to the hydraulic diffusivity with time going on when the system is approaching its steady-state condition. The flow system will become completely independent of the hydraulic diffusivity when it is at the steady-state condition. Therefore, the insensitivity to changes of hydraulic diffusivity (as the saturated thickness and specific yield will not change considerably), is an indicator that the flow system has passed its transiency period (in which the flow is sensitive to the hydraulic diffusivity), and is close to its steady state. In contrast to the findings for K x , small perturbations in the parameters / and S y (between 0.0001 and 0.0005) were found to change the water table height by several centimeters. The water table was most sensitive to recharge, which was varied as a percentage of the precipitation (10%, 20%, and 30%). These values produced drastic changes of up to meters in the water table height.
One initial PEST calibration (parameter estimation simulation) (Doherty 1994) was performed using one observation well from Texas because of the availability of seven observations during the simulation period. This initial calibration attempt was performed to obtain an idea of appropriate values for K, S y , and recharge. Then, the observation wells from Mexico and the U.S. were included to perform an improved second calibration (Supporting Information-07). There are five parameters that can be modified during calibration with observation heads: the three conductivity parameters (K x , K y , and K x ), specific yield (S y ) and recharge. Every time that a parameter was calibrated, the model was run again to proceed to consideration of the next calibration parameter. During this first PEST calibration, recharge was calibrated first, and then S y and K x were calibrated. The parameter that had the most influence over the calculated water table height was recharge, which was modified by a multiplier of 0.66 (calculated with PEST). While the method previously described to calculate recharge (CN) does not take into account the soil moisture conditions before recharge occurs, overestimation of recharge was expected; the PEST calibration accounts in this case for that excess of recharge as soil moisture and subtracts it from the initial uploaded recharge values. The next most significant parameter was the specific yield (S y ), which was assigned an initial value of 0.22 prior to calibration. After the PEST calibration, the value obtained for S y was 0.20, which yielded more accurate water levels when compared with the measurements from the observation wells.
Based on the results of the sensitivity analysis and the lack of response of the water table height to variations in K, the values adopted for K x , K y , and K x were the lowest values reported by Boghici (2002) and Castillo Aguiñaga (2000), which were described in the Methods section.
The calibration times selected were the time steps close to when observations were available (years 2006, 2008, 2011, and 2014). The heads reported for the 222 observation wells were compared to the calculated water levels for each time step (month), resulting in a 95% confidence interval of water level, a standard error of the estimate from 0.36 to 0.38 m, a root mean squared error ranging from 5.35 to 5.75 m, and a correlation coefficient of 0.99, which are good indicators of an efficient calibration process (Supporting information-08).

DISCUSSION
The first of the major findings obtained from modeling the APN aquifer is that stored water was depleted between 2000 and 2017. The second major finding is that prior to extensive development of the aquifer ("pre-2000"), the water convergence zone was probably located directly below the Rio Grande, whereas after development began to accelerate, this zone shifted under the influence of pumping from wells located near the riparian zones on both sides of the border. The third major finding is that this pumping proximal to the river captured recharge from the river. Today, most of the groundwater discharged into the Rio Grande comes from Mexico, in contrast to the pre-2000 period, when similar volumes of water were discharged into the river from both countries. Furthermore, in addition to the diminished baseflow to the river, pumping on the U.S. side derives a greater annual volume of captured recharge from the Rio Grande than pumping on the Mexico side.

Initial Conditions and Water Level Evolution
The pre-2000 water table was generated from observations made in September 1999 and January 2000 by the USEPA and TWDB (Supporting Information-07). These water levels were interpolated using Kriging for visualization (Figure 5a). The observations made in 1999 were assumed to represent the initial (pre-2000) conditions of the water table. The pre-2000 potentiometric surface generally follows the topographic slope, which is typical of unconfined aquifers. On the Mexican side, water flows from high terrain in Serrania del Burro and Lomerio Peyotes to the low-lying flatlands between El Moral, Piedras Negras, and Guerrero and later into the Rio Grande Valley flowing generally from the southwest to the northeast. On the U.S. side, groundwater flows from north to south down from the Anacacho Mountains into the Rio Grande Valley (Figure 5a).
According to the simulation results, the water table contours reveal a change in the water levels in comparison to the pre-2000 conditions. This shifting of the water table contour lines resembles a drop in the water table over time, as a consequence of the higher water extraction through pumping in the APN aquifer, which is not fully replenished by recharge. The spatially averaged cumulative drop in the water table for the entire aquifer was 0.76 m over 18 years, which is approximately 0.04 m/year. Additionally, the depth to the water table is greater where pumping has been substantial in the area between Piedras Negras and Villa Union on the Mexico side, and a consequent local depletion is observed, as shown in Figure 5b. The simulated drawdowns were highly variable spatially, with maximum drawdowns of 18 m during drought periods. There are some local areas where the water table recovered after cessation of pumping and water infiltration from streams. Among these is the area between the cities of Allende and Piedras Negras. The area located south of Piedras Negras suffered a strong local depletion of 14.8 m as a result of the large number of production wells operating in the area (Figure 6a).
The most rapid annualized drawdown rate, 0.06 m/ year, occurred in 2011, as a result of the drought (Figure 6b, average plot). In contrast, water levels were relatively stable in 2010, with a drawdown rate of only 0.015 m/year, as a result of the extraordinary rainfall events that occurred during April, June, and July 2010 (Figure 6b, average plot). Water table declines were compared between two regions of the APN aquifer that differed considerably in their pumping intensity. In a region with no pumping, located north of Eagle Pass, the water table declined by 0.4 m over 18 years. In contrast, in a region with intensive pumping, located south of Piedras Negras, the water table declined by 14.2 m over 18 years. Also noteworthy is the modeled exponential decay trend of the drawdown for the southern Piedras Negras region (Figure 6b density of pumping wells. This may be a result of the strong influence of the pumping.

Impact of Drought on Water Tables
The United States drought monitor was used to identify the worst droughts from 2000 to 2017 and evaluate the response of the water table to those droughts. The U.S. drought monitor is obtained from a compilation of several drought index indicators taken from different models, such as the standardized precipitation index and the Palmer drought severity index, among others. It is generated by the National Drought Mitigation Center, the U.S. Department of Agriculture and the National Oceanic and Atmospheric Association (Svoboda et al. 2002). Detailed drought monitor outputs for Texas were used to identify the driest months for Southeast Texas: May 2006, May 2009, September 2011, and April 2013 According to the drought monitor, Maverick and Kinney Counties were not under drought in December 2000. This was also one of the periods of least pumping activity in the APN aquifer. Water table depths were compared for the periods of December 2000 and September 2011 to evaluate the effects of exceptional droughts on the water table across the APN aquifer (Figure 7).
Most of the water table depth changes were observed at a local scale, where excessive pumping has caused the water table to fall between the regions of Piedras Negras and Guerrero, reaching depths of close to 18 m. The depletion seen in the region between Allende and La Union, however, is probably a consequence of minimal recharge, rather than excess water extraction, because few wells are located in the area.

Water Flow across the Rio Grande System
According to the simulation results, the groundwater flow directions in the APN aquifer follow the land topography and converge in the Rio Grande. The groundwater in the Texas side of the aquifer flows from the Anacacho Mountains to the southwest, and the groundwater in Coahuila flows from the highlands of Serrania del Burro and Lomerio Peyotes to the northeast, toward the Rio Grande ( Figure 8a). Under early-2000 conditions, the groundwater convergence zone would have been located underneath the Rio Grande, as explained by Williams (1993), where the groundwater discharges to the river and vice versa, creating a water mixing zone. Pumping may have altered the timing, volumes, and direction of fluxes between the river and the aquifer.
A detailed evaluation of the modeling results for early-2000 conditions in the Rio Grande riparian zones and groundwater fluxes into the Rio Grande indicates that the model suggests that the river is likely a "gaining river," specifically in the regions to the north of Quemado, Texas and El Moral, Coahuila (Figures 8b, 8c, and 9, cross section A-A 0 ). In contrast, the river becomes a "losing river" in the southern region near Piedras Negras, Coahuila, where the river recharges the aquifer during dry seasons, as seen in Figures 8b, 8c, and 9 (Cross section B-B 0 ). According to the model, however, the river is predominantly a gaining river because the water table contour lines are pointing upstream along most of the course of the Rio Grande (Winter 1998) (Figure 8c). A local effect of the pumping wells on river flows is apparent (Figure 9, cross sections C-C 0 and D-D 0 ). The groundwater converges on one side of the river but not below, as it does under early-2000 conditions. Depending on the proximity of the pumping wells to the river and their pumping rates, the groundwater convergence zone will shift toward the region with greater pumping rates. In Figure 9 (cross section C-C 0 ), the groundwater flow switches to the Texas side because of pumping from the aquifer in South Quemado. In some cases, groundwater from the APN aquifer discharges into the Rio Grande, and afterward, the water returns to the APN aquifer downstream. This phenomenon of water discharge and later infiltration back into the aquifer was referred to as flow-through by Hoehn (1998).
In the area located north of Piedras Negras, Coahuila, the density of pumping wells in Mexico is higher than on the Texas side (Figure 8b). This causes the groundwater convergence zone to move into Mexico and causes flow-through from the east to the west, as seen in Figure 9 (cross section D-D 0 ). The main implication of this alteration in the water flow patterns is that this type of alluvial aquifer, wherein the river acts as a geopolitical border, has "little transboundary flow" (Rivera 2015). The exception to this condition is if the extraction of groundwater affects the baseflow of the river enough to induce substantial changes in hydraulic heads, as observed in this study. This transforms the aquifer system into a transboundary groundwater flow system (Rivera 2015).

Water Budget
A water budget was generated for the areas described previously in the subsection titled "Water budget zones," with the aim of determining the groundwater discharge into and recharge from the Rio Grande to the separate parts of the aquifer located in Mexico and the U.S. The estimated water budget for early-2000 conditions is shown in Figure 10.
The water budget appears unbalanced because there are still unknown pieces, such as the groundwater flow entering the aquifer and the change in storage. The APN aquifer recharge is dominated by infiltration of rainfall, but because of the limited number of rainfall events and the arid to semiarid climate, substantial amounts of water are lost from the system through ET. Furthermore, even when water exchange between the Rio Grande and the riparian zones were observed in the water budget analysis, the river was an important discharge component of the APN aquifer. Figure 11 shows the total amounts of water flowing into and out of the APN aquifer over the 18-year simulation. The data were obtained from the accumulated mass balance calculation in VMF. The mass balance indicates a net loss of water, meaning that the outflow was higher than inflow in the APN aquifer, depleting groundwater stored in the aquifer. The inflow sources, such as recharge from rainfall and streams, are not sufficient to meet the demand from the aquifer, and the system has been depleted to some degree.
Theoretically, in an aquifer under pre-development conditions, the inflows and outflows should be equal. Under development conditions, however, activities such as water pumping, surface water capture, and land use changes affect the amounts of water entering and leaving the system. Under these dynamic conditions, storage becomes a source of water (Alley et al. 1999). It is common to observe a change in storage at the beginning of pumping in any aquifer because the system responds to the water withdrawal. Water flowing from storage will cease when the system reaches a new equilibrium and the inflows equal the outflows again, provided that new sources of recharge, such as a perennial river, can be accessed by the aquifer in response to lowering of the water table (Alley et al. 1999).
A comparison of the annual inflows and outflows to and from the APN aquifer showed that storage amounts from inflow and outflow differed throughout the simulation period. The implication of this discrepancy is that the aquifer has not reached equilibrium and instead that part of the water moving out of the aquifer through the rivers, ET, and wells was derived from storage. Another important observation is that most of the water extracted, 95%, was extracted by the Mexican wells; only 5% was extracted by wells located on the U.S. side. According to Figure 12, in most of the years during the simulation period, the outflows exceeded inflows, with the greatest negative difference of À10.9% occurring in 2008. There were years in which the differences between inflows and outflows were slightly positive (2004, 2005, and 2010), ranging from 0.03% to 0.41%.
The Rio Grande is an important source of water infiltration into the APN aquifer. As discussed in previous sections, areas located in the southern region of the Rio Grande function mainly as recharge zones, whereas those in the northern region function mainly as discharge zones. However, the Rio Grande is a net discharge zone for the APN aquifer. Under pre-2000 conditions, the river received similar volumes of groundwater from the U.S. and Mexico portions of the aquifer. Under the influence of pumping wells in the riparian zones, however, groundwater discharge into the Rio Grande has been modified. According to Figure 13, the cumulative percentages of water volumes in December 2017 were distributed as follows: most of the water discharged into the Rio Grande came from Mexico (55.5%), whereas the remaining 44.5% came from the U.S.; under postdevelopment conditions, 51.1% of the water infiltrated from the Rio Grande recharged the U.S. portion of the APN aquifer, and the remaining 48.9% recharged the Mexico portion.

Limitations and Future Work
A limitation of this study is the lack of information that was available on the U.S. portion of the aquifer, which prevented modeling of the initial conditions in this region at a high spatial resolution.
Ideally, to determine infiltration through streams precisely, model calibration should include observed heads and flows. Although calibration was performed using information from observation wells, it was not possible to calibrate the model using flows because of the lack of this type of data. Therefore, caution should be used in applying this model to seepage analysis.
The main reason for averaging ET and recharge values by assigning monthly averages to the entire area was to simplify the model and thereby reduce running and calibration times and the possibility of nonconvergent solutions. This averaging reduces the spatial resolution of the datasets. In future work, we intend to integrate detailed ET and recharge information and thereby avoid averaging over large areas.
On part of the future work to be conducted will focus on developing a more detailed hydrogeological model. This will be achieved by augmenting several of the datasets developed from remote sensor data using field data with higher temporal and spatial resolutions, such as rain gage data instead of satellite images. Refining the grid in the APN aquifer, specifically in the transect between El Moral-Quemado and Piedras Negras-Eagle Pass, would allow a more detailed analysis of the pathways for water flow near the geopolitical boundary between Mexico and the U.S. Running the numerical model with cross-formational flow taken into account would refine flowpaths, determine other sources of recharge rather than precipitation, and provide more accurate water budget calculations. In addition, a detailed water budget  could be generated from this refined numerical model that would more accurately calculate the water flows moving between the aquifer zones, as well as the interactions between the Rio Grande and the surrounding zones of the aquifer. A more accurate calculation of recharge, taking into account land use, would also improve the model's ability to predict the effects of future pumping activities and land use changes on the behavior of the river-aquifer system. However, while VMF is not an ideal tool for integrating surface water phenomena with groundwater flow, the variation in recharge related to land use can be modeled using specialized watershed software, and the resulting recharge can be included in the groundwater model as recharge.
The generation and analysis of hydrographs and comparison with hydrogeochemical data will further refine the findings obtained concerning surface water-groundwater interactions between the Rio Grande and the APN aquifer.

CONCLUSIONS
The following conclusions were drawn from the results of this study: 1. A spatially averaged cumulative drawdown of 0.76 m was calculated over the 18-year simulation period. Even though the average drawdown rate lessened after 2011, the aquifer has not recovered to pre-2000 water levels. 2. Under pre-2000 conditions, the water flowpaths from the aquifer converged into the Rio Grande, but after the number of pumping wells installed adjacent to the river began to increase considerably, the flow convergence zone shifted to the Mexican side of the border, because of greater combined pumping rates. This modification of the baseflow of the river and the change in hydraulic heads allows us to classify the APN aquifer as a transboundary groundwater flow system (Rivera 2015). 3. According to the water budget, the amounts of water extracted from the aquifer surpassed the inflows, indicating that stored water is being depleted. This finding indicates that the APN aquifer is over-exploited, because water extracted naturally or artificially is not entirely recovered on a long-term basis. 4. Mexico pumps 95% of the total groundwater volume extracted per month from the APN, whereas Texas only pumps 5%. This reflects the dependency of Mexico on the APN aquifer. The inhabitants and industries located on the Texas side of the border rely mostly on the deeper Edwards aquifer. 5. The portion of the aquifer on the U.S. side of the border is discharging less water (44%) than the portion on the Mexico side (56%). However, the proportion of recharge to the APN aquifer from the Rio Grande is higher on the U.S. side (51%) than on the Mexico side (49%).

SUPPORTING INFORMATION
Additional supporting information may be found online under the Supporting Information tab for this article: Detailed description of geology, maps, and tables produced during data analysis, conductance calculation, calibration outputs, and hydrochemical maps.