**2. Materials and Methods**

### *2.1. Study Area*

The Los Laureles Canyon Watershed (LLCW) is a transboundary urbanizing catchment located in the northwestern part of Tijuana, Mexico, which flows into the Tijuana River Estuarine Reserve, USA (Figure 1). The total catchment area is 11.6 km2, with 93% in Mexico and 7% in the USA. The climate in LLCW is Mediterranean, with a rainy season during the winter and annual mean rainfall of approximately 240 mm. Most of the erosional storm events occur during the winter. The regional geology (San Diego formation) includes marine and fluvial sediment deposits of conglomerate, sandy conglomerate, and siltstone. Soils are sandy with a wide range of cobble fraction, and are dominated by abrupt slopes (15◦, average), which encourages gully formation and results in high erosion rates.

The LLCW is an uplifted and incised marine terrace, where the soil types are controlled in part by the underlying geology. The conglomerate geology in the northern part of the watershed has steep, competent valley walls with relatively flat buttes and mesas. In the central part of the watershed, the sandy conglomerate geology with low cobble fraction has lower slopes and rounded hilltops. The southern part of the watershed is a relatively flat, non-incised conglomerate. A narrow valley floor has Quaternary alluvium, but most of this has been paved or channelized.

Land use in the LLCW is predominantly mixed urban and rural, and was urbanized starting in 2002 with many illegal housing developments ("invasiones"). Gullies form on unpaved roads, which affect civil infrastructure in the upper watershed [31] and downstream ecosystems [32]. Such gullies are filled in with sediment following storms, and this management practice should be taken into account in developing the sediment budget and in soil erosion modeling for the watershed.

Sediment from the LLCW has buried native vegetation in the Tijuana River Estuary, which is located downstream of LLCW in the United States. In response, sedimentation basins were constructed at the outlet in the US in 2004, which costs \$3 million USD to clean annually [32]. Stormflow and erosion also threaten human life, which causes damage to roads and houses in Tijuana [31]. The primary sources of sediment from LLCW are gully formation on unpaved roads, channel erosion, and sheet and rill erosion from unoccupied lots in Tijuana [8].

#### *2.2. Field Data Collection and Model Setup*

A summary of the data collection activities is reported in Reference [33]. Briefly, a tipping-bucket rain gauge station (RG.HM in Figure 1) was installed in February 2013. A pressure transducer (PT) (Solinst, water level logger) was installed in a concrete channel at the watershed outlet in December 2013 to record the water stage at 5-minute intervals (Figure 1). The stage-discharge relationship was determined using Manning's equation and flow velocity measurements. Manning's roughness coefficient (*n*) was based on field measurements of discharge in 2016 and 2017, which was used to back-calculate a Manning's *n*. The discharge measurements were also used to create a stage-discharge relationship for a stream gauge in the US (RG.GC) to complete our observations when the pressure transducer malfunctioned. The PT data were also validated and supplemented using time-lapse photographs of the water stage at the PT station. Suspended Sediment Concentration (SSC) measurements were taken at 10 different locations during a storm on 27 February, 2017 to explore spatial patterns of sediment production within the watershed. Annual sediment load data was collected in

two large sediment traps at the watershed outlet. Data on the quantity of sediment removed annually (2006–2012) from the traps were available from the Tijuana River National Estuarine Research Reserve (TRNERR), corrected for trap efficiency, and used for model calibration (Figure 1).

**Figure 1.** Location of the Los Laureles Canyon Watershed, main channels (SW, Main, and SE), and monitoring stations including sediment traps and rain gauges (RG). Inset shows the geographic locations of nearby rain gauges used in this analysis to span the rainfall time series.

A map of soils and associated parameters needed for the model were not available for the watershed, so soil type and parameters were mapped by modifying an existing geology map [34] and a correlation between geology and soil type that was established using the Soil Survey Geographic Database (SSURGO) from the United States Department of Agriculture. The geology map was modified based on field data collected during September 2015 and visual interpretation of high resolution imagery on Google Earth. First, a seamless cross-border geological map was created using the Instituto Metropolitano de Planeación de Tijuana, Baja California Mexico (IMPLAN) [34] geology map for Mexico and a geology map from the US [35]. The US soils that occurred on geologic types found in LLCW were identified as candidate soils for the study watershed. See Biggs et al. [33] for a full description of field and laboratory data collection.

Three main soil types were identified: (1) Los Flores formation (Lf) is a loamy fine sand, (2) the Chesterton formation (CfB) is conglomerate dominant with a fine sandy loam matrix, and (3) the Carlsbad formation (CbB) is a gravelly loamy sand. Once candidate soils were determined from the US geology and soils maps, the SSURGO soil characteristics were extracted for all horizons for comparison with data collected in the field. Soil samples were then collected from different geology types and analyzed for texture to compare with the SSURGO database and with the observed soil texture in the sediment traps at the watershed outlet (Figure 2). Samples (*N* = 25) were collected from road cuts and other exposed profiles from the near-surface (10–50 cm) and from the subsurface (>50–100 cm). The cobble percentage was determined through point counts along a 1 m transect through each distinct horizon. A bulk sample of sediment smaller than coarse gravel (<32 mm) was collected for texture analysis, analyzed in the laboratory using dry sieving to separate a 2-mm fraction and the pipette method for fines (<2 mm). Soil texture for all soil samples collected in LLCW and near the US-Mexico border was plotted in ternary diagrams and compared to SSURGO surface and subsurface soil texture (Figure 2). For each soil group, the SSURGO soil types that most closely matched the mean texture from the soil samples, were selected and used to update the soils map for LLCW. For some areas, the texture from the samples was similar to SSURGO data (northern part of the watershed). For the CfB soils, the texture from the samples did not match the SSURGO texture (southern part of the watershed), so a new geologic type (CfB.MX) was created, supported by field and laboratory data, to include in the AnnAGNPS model. Lastly, polygons delineating soil types were created by first determining the relationship between soil color, landform, and soil type for soils in the US, and then extrapolating those relationships to map similar soils in the LLCW (Figure 3).

**Figure 2.** Ternary diagram of the mean grain sizes observed in the surface (SURF) and subsurface (SUB > 50 cm depth) soil layers for each soil type in the watershed, and in the sediment traps in Mexico (MXSB) and at the outlet in USA (USTRAP).

**Figure 3.** Geology map (**a**) and updated soils map for LLCW (**b**). Las Flores (Lf), fine sandy loam, dominates the central portion of the watershed (orange). CfB.MX represents the Chesterton sandy loam (CfB), but with a cobbly surface horizon. Carlsbad (CbB) and "CfB.US" soils extend south from the US/Mexico border.

The critical shear stress (τ*c*) and soil erodibility in the non-cobbly sandy conglomerate (Lf) were taken from Gudino-Elizondo et al. [14], who used a mini-jet erosion test device following the methodology described by Hanson [36]. Values of τ*<sup>c</sup>* for conglomerate soil were initially taken from USGS [37] and were modified during calibration.

The AnnAGNPS model requires daily precipitation data. To extend the simulation period to include the period after the installation of the sediment trap (2004–2017), the precipitation data collected from February 2013 to 2017 at RG.HM were compared with rainfall data from nearby stations in the United States (Figure 1) to select the best rain gauge to use for rainfall data from 2004–2013.

Application of AnnAGNPS can be challenging in a watershed with steep topography and sediment coarser than sand (Figure 4) because the model does not simulate mass wasting processes, only transports sediment up to coarse sand (2 mm), and is designed for mixed-use watersheds in agricultural areas. Mass wasting, including shallow landslides, was observed in the study area, and coarse sediment accounts for ~10–15% of the sediment in the traps at the watershed outlet (unpublished data). However, the valley floor at the base of the steep slopes that are most likely to experience landslides, has been graded and paved for roads on either side of the channel, which limit the transport of coarse material to the channel from landslides or other hillslope processes. Field observations suggest that landslides typically terminate on these flat road segments or other graded areas, and the coarse material that accumulates at the toe of a landslide is periodically cleared mechanically. In this scenario, we assumed that all coarse material is from the channel. The sediment load from channel erosion was taken from Reference [38] and added to the total modelled hillslope erosion. Lastly, the RUSLE equation to estimate sheet and rill erosion was designed for relatively flat agricultural hillslopes, and may not be valid for steep hillslopes (>30%) [39]. In this case, we assumed that the model application was valid in most of the study watershed because 87% of the total watershed area has a slope gradient of less than 30% (mean slope = 15%).

The spatial variability of topography, land cover, soils, and management properties within the catchment area was represented in the model by discretizing the watershed into cells that are relatively similar in slope, soils, and land use. A LIDAR-derived Digital Elevation Model (DEM) (3 m, horizontal resolution) of the LLCW (sponsored by the County of San Diego California, USA) was used as input for a topography-based method (TopAGNPS) [26]. The TopAGNPS method was used to (i) delineate surface flow paths, (ii) subdivide the total catchment area into sub-catchments (cells) along drainage segmentations, and (iii) estimate representative cell parameters, such as slope, area, and soil and management attributes. Cells sizes were based on user-defined values of the Critical Source Area (CSA), which is the minimum drainage area required to form a channel, and the Minimum Source Channel Length (MSCL), which prunes the channel network of channels shorter than the specified MSCL value.

To characterize the hillslope and reach units within the AnnAGNPS model, a CSA of 1 ha and a MSCL of 50 m were assigned, based on field observations, to characterize the hill-slope and reach units within the AnnAGNPS model. LLCW was discretized into 1147 sub-catchments (AnnAGNPS cells) and 462 channels (AnnAGNPS reaches). The cell sizes ranged from 9×10−<sup>6</sup> to 0.1 km2 (Figure 4).

TopAGNPS was also used to map Potential Ephemeral Gullies (PEGs) throughout the LLCW following the methodology described by Momm et al. [26]. This method provides an automated estimate of the downstream-most locations of knickpoints (i.e., PEGs), which are used within AnnAGNPS to calculate the length of ephemeral gullies in the landscape. The approach uses improvements on the EGEM described in Gordon et al. [40] and, more recently, revised by Bingner et al. [12]. In the model, the gullies are filled in once per year at the end of each wet season, which corresponds to observed management practices, but may under-represent the filling frequency on larger main roads that are typically filled between each storm event that generates gullies.

The watershed hydrology module of AnnAGNPS uses the SCS Curve Number method [27] to estimate storm event runoff from precipitation in each cell. The storm event water peak discharge and the time-to-peak are determined for the hydrograph at each reach section and at the watershed outlet, following the TR-55 [27] approach that utilizes the time of concentration for each cell and reach, determined with TOPAGNPS, total daily runoff determined from AnnAGNPS, and the storm type entered as an input parameter [12].

**Figure 4.** (**a**) AnnAGNPS cells, reaches, mainstream channels, and (**b**) slope gradient at the Los Laureles Canyon watershed (LLCW).

The spatial resolution of the DEM has little impact on AnnAGNPS runoff volume estimates [11,41], but soil erosion and sediment loads can change with DEM resolution, since resolution impacts slope [42]. LIDAR-derived 3-m DEMs should improve the model performance in the study watershed compared with applications that use more-commonly available 10 m or 30 m resolution DEMs. Initial AnnAGNPS reaches did not follow the road network, so the road segments were "burned" into the DEM by lowering the elevation on the DEM cells falling on roads by 1 m.

A land use map was created using Google Earth (11 November 2012, 2017 Digital Globe) imagery based on visual classification of seven categories (rangeland, highway road, paved residential roads, dispersed urban unpaved, unpaved rural roads, unpaved residential roads, and sediment basin). The accuracy of the land use data was validated by comparing land use categories with ground-based photography and field data collection. The land use map was overlain on the AnnAGNPS cells to populate the required hydrologic and management parameters needed for AnnAGNPS. The soils map was used to link the required physical variables from the SSURGO database to the model such as soil texture and erodibility, bulk density, and saturated conductivity. Tillage depth is the depth to an impervious soil layer, which limits the potential depth of gullies, and was determined as the depth of the gullies observed in the field [14]. The main equations solved within AnnAGNPS to estimate soil erosion are listed in Table 1. A detailed description of these equations are given by Bingner et al. [12].


Gully width W = 9.0057 ×

Head-cut migration erodibility coefficient Kh <sup>=</sup> 0.0000002/ <sup>√</sup>

Annotation: Sy = sediment yield by sheet and rill erosion (Mg/ha): Q = surface runoff volume (mm). qp = peak rate of surface runoff (mm/s) and K, L, S, C, P are RUSLE factors. W = gully width (m): Qp = peak discharge at the gully head (m3/s). S = the average bed slope above the gully head (m/m). Kh = head-cut migration erodibility coefficient (m3/s/N): τ*<sup>c</sup>* = the critical shear stress (N/m2).

See Gudino-Elizondo et al. [14] for a detailed description of the usage of these equations in the study watershed.
