Simulation of Nitrogen Migration and Output Loads Under Field-Scale in Small Watershed, China

1. Introduction

Small watersheds are the fundamental unit for watershed water resources management and play an irreplaceable role in controlling soil erosion, managing pollution sources, and improving water environments [1]. Under the profound influence of agricultural intensification and land-use changes, small watersheds have become key areas for controlling agricultural non-point source pollution [2,3]. According to the results of the second national pollution source census, nitrogen emissions from agricultural sources account for 46.5% of the total emissions in China. The migration and transformation processes of nitrogen have a decisive impact on the eutrophication of receiving water bodies [4,5]. Thus, there is an urgent need to conduct detailed studies at the source scale.

Although significant progress has been made in nitrogen migration research at the watershed scale, physical process models such as SWAT and HSPF have shown good performance in large-scale applications [6]. However, these models have limitations in precisely capturing the details of nitrogen migration at the field scale and quantifying the associated output loads [7]. Nitrogen migration at the field scale exhibits high spatial heterogeneity and is strongly influenced by microtopography (such as field embankments and ditches) and intricate land-use patterns. Thus, models capable of adapting to small-scale complex topography are needed. The SCS-CN model, due to its minimal parameter requirements and computational simplicity, has become an effective tool for rainfall-runoff simulation, particularly in agricultural areas with complex land-use types [8,9]. The MUSLE model, by incorporating runoff factors, improves the traditional USLE model, enabling more accurate estimates of soil erosion and pollutant migration [10,11].

High-precision Digital Elevation Models (DEMs) are critical for simulating nitrogen migration at the field scale, as their resolution directly affects the accuracy of microtopography feature identification and hydrological process simulations. The Multi-Flow Direction (MFD) algorithm, which overcomes the limitations of traditional single-flow direction algorithms, demonstrates significant advantages in complex topographies, especially when combined with high-precision DEMs [12,13,14]. It significantly enhances the accuracy of hydrological process simulations.

Field-scale nitrogen migration is influenced by various factors, including soil type, crop configuration, microtopography, and management practices [15]. Research has shown that different agricultural planting patterns significantly alter nitrogen loss pathways and efficiency [16,17], while vegetation coverage and precipitation characteristics further regulate nitrogen migration intensity in small watersheds [18]. However, existing research primarily focuses on larger scales, and there is still insufficient quantitative analysis of nitrogen migration mechanisms at the field scale [19,20]. Therefore, further exploration of these mechanisms is necessary [21,22].

The Jianghuai watershed region is an important agricultural production base in China, facing dual ecological pressures from the Yangtze and Huaihe River basins. Non-point source pollution is a significant issue in this region, which features complex topography and diverse land-use types, making it an ideal location to study nitrogen migration mechanisms at the field scale. This study focuses on a typical agricultural small watershed in the Jianghuai watershed. We use UAV photogrammetry technology to construct a centimeter-level spatial database, coupling the SCS-CN and MUSLE models to quantify nitrogen output loads, and apply the MFD algorithm to simulate nitrogen migration trajectories. By analyzing the synergistic effects of topography and land use on nitrogen migration, this study provides a basis for the fine-scale control of non-point source pollution in small watersheds.

2. Materials and Methods

2.1. Study Area Overview

The Jianghuai Watershed region (31°25′–32°40′ N, 116°50′–118°10′ E) is located at the boundary between the Yangtze River and Huaihe River basins, holding significant geographical and ecological importance. The region features a notable topography, with typical microtopographic features (such as field embankments and ditches), making it an ideal area to evaluate the role of high-resolution DEMs in simulating nitrogen migration at the field scale. This study focuses on a typical small watershed in an agricultural irrigation area within the Jianghuai Watershed (Figure 1), with a watershed area of approximately 1.27 km², characterized by diverse land-use types and complex nitrogen sources. The small watershed exhibits typical characteristics of agricultural non-point source pollution in the Jianghuai Watershed, including greenhouse planting, traditional farmland, and forest land, providing a good sample for studying nitrogen migration patterns under different land-use types. In recent years, due to agricultural intensification and other human activities, land-use types and agricultural structures in this region have undergone significant changes, leading to an intensification of non-point source pollution. Despite the implementation of various pollution control measures, the non-point source pollution situation remains severe and requires urgent management.

2.2. Technical Approach

This study selected this region to reveal the spatial distribution patterns of nitrogen migration at the field scale through the use of high-precision DEMs and multi-model coupling methods, aiming to provide technical support for fine-scale pollution control in similar agricultural small watersheds. The technical approach is outlined in Figure 2.

The study is divided into four main components: First, UAV photogrammetry, field sampling, and experimental analysis are employed to obtain centimeter-level high-resolution data and chemical indicator data; the data is preprocessed to generate orthophotos, land-use data, DEMs, and nitrogen data for soil and water bodies; an SCS-CN model is coupled with the MUSLE model to quantify nitrogen output loads based on nitrogen concentration; nitrogen migration trajectories are simulated to reveal the paths of nitrogen migration in the watershed; and the multi-flow algorithm is used to simulate the migration process of different nitrogen forms (TN, NH₄⁺-N, NO₃⁻-N). The primary data sources for this study are listed in Table S1.

2.3. Research Methods

2.3.1. Construction of a Centimeter-Level Database

• Acquisition of High-Precision DEM:

In this study, data collection was conducted using the DJI Phantom RTK drone. The photogrammetry method employed is based on existing literature, including flight plan design, drone photogrammetry experiment design (with a flight height of 50 m, flight speed of 5 m/s, longitudinal overlap set at 80%, and lateral overlap at 70%), and control measurement procedures. The images collected were processed using DJI Terra v2.4.1 software to generate orthophotos and DEMs, with a spatial resolution of 0.1 m.

2.

Land Use Interpretation

Based on the high-precision orthophotos of the study area, ArcGIS 10.8 software (Esri, USA) was used in conjunction with field surveys and visual interpretation to classify the land-use types in the study area. The land-use types were divided into 14 categories, resulting in a high-resolution land-use map (0.1 m × 0.1 m). The specific land-use types are shown in Figure 3(b).

3.

Collection and Analysis of Watershed Nitrogen Concentration Data

In August 2022, systematic field sampling was conducted in the study area. A total of 15 water quality monitoring points and 49 soil sampling points were set up. Water and soil samples were delivered to the laboratory within 24 hours after collection. In the laboratory, all samples were stored in professional refrigerated equipment at 4°C until the nitrogen concentration measurements were completed. Total nitrogen concentrations were measured using the alkaline potassium persulfate digestion UV spectrophotometric method (HJ636-2012); nitrate nitrogen concentrations were measured using the UV spectrophotometric method (GB/T 5750.5-2023); ammonium nitrogen concentrations were measured using the salicylic acid-hypochlorite spectrophotometric method (GB/T 7481-1987). The results were processed using ArcGIS 10.8 software to perform spatial interpolation of soil and water samples, constructing a spatial database for nitrogen concentrations in the watershed.

2.3.2. SCS-CN Model Runoff Calculation

The SCS-CN runoff model was used in this study to quantitatively calculate the runoff of the watershed. The model is based on the principle of water balance and introduces the curve number (CN) as a key parameter, which integrates the impact of various land-use types on the precipitation-runoff conversion process in the watershed [23].

The core formula is:

$$\left\{\begin{array}{c}Q={\displaystyle \frac{{\left(P-I_a\right)}^2}{P-I_a+S}};P>I_a\\ Q=0;P\le I_a\end{array}\right.$$

(1)

Where:

$$I_a=\mathsf{\lambda }\mathrm{S}$$

(2)

$$\mathrm{S}=25400/\mathrm{CN}\text{–}254$$

(3)

Where: Q is the runoff (mm); P is the precipitation (mm); Ia is the initial loss (mm).

λ is the initial loss coefficient, typically taken as 0.2; S is the potential maximum retention (mm); CN is the curve number (dimensionless, ranging from 0 to 100). In the calculation, CN values for different land-use types in the study area were obtained from related literature in Anhui Province and surrounding areas. Daily precipitation data were used to simulate the daily runoff process in the watershed, and the total annual runoff was then calculated [24,25,26].

2.3.3 MUSLE Model Sediment Yield Calculation

The MUSLE (Modified Universal Soil Loss Equation) model is an improved soil erosion prediction model based on the traditional USLE model, developed by the U.S. Department of Agriculture. The model improves prediction accuracy for erosion from single rainfall events by introducing the runoff factor to replace the original rainfall erosivity factor, making it a vital tool for soil erosion research in small watersheds [27].

The basic equation for the MUSLE model is:

$$\mathrm{A}=11.8\text{×}{\left(\mathrm{Q}\text{﹒}{q}_p\right)}^{0.56}\text{×}\mathrm{K}\text{×}\mathrm{LS}\text{×}\mathrm{C}\text{×}\mathrm{P}$$

(4)

Where: A represents the soil intensity from a single rainfall event (t/km²); Q is the total runoff (mm); q_p is the peak flow (m³/s); K is the soil erodibility factor; LS is the topographic factor (where L is the slope length factor and S is the slope factor); C is the vegetation cover and management factor (dimensionless, ranging from 0 to 1); P is the conservation practice factor (dimensionless, ranging from 0 to 1).

1.

Runoff Volume Q and Peak Flow q_p

In the formula,

$$q_p=C\times i\times F$$

(5)

$$C=Q/P$$

(6)

Where: C is the runoff coefficient; i is the rainfall intensity; F is the watershed area.

Peak flow values for different land-use types were calculated using daily precipitation data obtained from the European Center for Medium-Range Weather Forecasts and formulas (5) and (6), and assigned values using ArcGIS software [28].

2.

Soil Erodibility Factor K

Soil erodibility reflects the sensitivity of soil to erosion. The K value is calculated using the equation from the EPIC (Erosion/Productivity Impact Calculator) model, as proposed by Williams [29]:

$$\begin{array}{cc}K=& \{0.2+0.3\mathit{exp}[-0.0256SAN(1-{\displaystyle \frac{SIL}{100}}]\}\times \\ & {\left({\displaystyle \frac{SIL}{CLA+SIL}}\right)}^{0.3}\times (1.0-{\displaystyle \frac{0.25C}{C+exp(3.72-2.95C)}})\times \\ & (1.0-{\displaystyle \frac{0.7SN}{SN+exp(22.9SN-5.51)}})\end{array}$$

(7)

Where: CLA, SIL, and SAN are the percentages of clay, silt, and sand in the soil; C is the organic carbon content in the soil; SN = 1 − SAN/100.

The K values for the different land-use types in the study area were determined by referring to the China Soil Database and other relevant studies, and the results were mapped to the land-use type map.

3.

Slope Length and Slope Factor LS

Slope and slope length are topographic factors that affect soil erosion. The longer and steeper the slope, the higher the soil erosion rate. The slope length factor L was calculated using the formula proposed by Desmet et al. [28,30]:

$$L={\left({\displaystyle \frac{\lambda }{22.1n}}\right)}^n,n=\left\{\begin{array}{c}0.2\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}(\theta ⩽1^{\circ })\\ 0.3\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}(1^{\circ }<\theta ⩽3^{\circ })\\ 0.4\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}(3^{\circ }<\theta ⩽5^{\circ })\\ 0.5\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}(\theta >5^{\circ })\end{array}\right.$$

(8)

$$\lambda =Flowaccumulation\times Cellresolution$$

(9)

Where: λ is the flow accumulation and cell resolution; n is the slope length index.

The slope factor S is calculated as:

$$S=\left\{\begin{array}{cc}10.8\mathrm{s}\mathrm{i}\mathrm{n}\theta +0.03& \theta <5^{\circ }\\ 16.8\mathrm{s}\mathrm{i}\mathrm{n}\theta -0.5& 5^{\circ }<\theta \le {10}^{\circ }\\ 20.204\mathrm{s}\mathrm{i}\mathrm{n}\theta -1.2404& {10}^{\circ }<\theta \le {25}^{\circ }\\ 29.585\mathrm{s}\mathrm{i}\mathrm{n}\theta -5.6079& \theta >{25}^{\circ }\end{array}\right.$$

(10)

Where: θ is the slope (°); λ and α are parameters related to the slope length.

Slope length and slope values were extracted from the 0.1m DEM data in ArcGIS and used to calculate the topographic factor LS.

1.

Vegetation Cover and Management Factor C

The vegetation cover and management factor C is influenced by vegetation growth, type, and coverage. This study assigned values to different land-use types based on research from surrounding watersheds [30].

Table 1. Vegetation Cover and Management Factor C Values for Different Land Use Types.

Land-Use Type	Water Body	Forest Land	Shrub Land	Residential Area	Bamboo Forest	Grassland	Nursery Land	Vegetable Field	Bare Land	Greenhouse	Dryland	Paddy Field	Road
C value	0.00	0.01	0.02	0.02	0.02	0.06	0.08	0.09	0.15	0.25	0.25	0.25	0.40

Table 1. Vegetation Cover and Management Factor C Values for Different Land Use Types.

Land-Use Type	Water Body	Forest Land	Shrub Land	Residential Area	Bamboo Forest	Grassland	Nursery Land	Vegetable Field	Bare Land	Greenhouse	Dryland	Paddy Field	Road
C value	0.00	0.01	0.02	0.02	0.02	0.06	0.08	0.09	0.15	0.25	0.25	0.25	0.40

2.

Conservation Practice Factor P

The P value, which indicates the level of soil conservation, was assigned based on the land-use situation in the study area and surrounding areas [31,32].

Table 2. Conservation Practice Factor P for Different Land Use Types.

Land-Use Type	Water Body	Paddy Field	Grassland	Vegetable Field	Bare Land	Greenhouse	Dryland	Road	Forest Land	Bamboo Forest	Shrub Land	Nursery Land	Residential Area
P value	0.00	0.01	0.20	0.20	0.20	0.20	0.25	0.50	0.50	0.50	0.60	0.60	1.00

Table 2. Conservation Practice Factor P for Different Land Use Types.

Land-Use Type	Water Body	Paddy Field	Grassland	Vegetable Field	Bare Land	Greenhouse	Dryland	Road	Forest Land	Bamboo Forest	Shrub Land	Nursery Land	Residential Area
P value	0.00	0.01	0.20	0.20	0.20	0.20	0.25	0.50	0.50	0.50	0.60	0.60	1.00

2.3.4. Nitrogen Migration Direction and Amount Calculation Based on Multi-Flow Direction Algorithm

The Multi-Flow Direction (MFD) algorithm is used for flow distribution in hydrological models to address the limitations of traditional single-flow direction algorithms (such as D8), especially in areas with flat terrain or gentle slopes. The D8 algorithm limits the flow of material to the steepest slope, which may cause unrealistic distribution in flat or gently sloped regions. In contrast, the MFD algorithm allows material to flow to multiple neighboring pixels, providing a more accurate simulation of material distribution in complex topographies [33].

This method is particularly suitable for processing large-scale Digital Elevation Models (DEMs), as it better simulates water flow across the entire surface, especially in areas where flow disperses and converges. The slope s_ij from pixel i to pixel j is calculated as:

$$s_{ij}={\displaystyle \frac{dem\left(i\right)-dem\left(j\right)}{d_{ij}}}$$

(11)

Where: s_ij is the slope between pixels i and j; dem(i) and dem(j) are the elevations of pixels i and j; d_ij is the distance between pixels (usually the unit grid size).

Next, flow A_i is distributed to neighboring pixels based on the slope:

$$A_i=\sum _{j\in N\left(i\right)} {\displaystyle \frac{A_j.S_{ij}}{\sum _{k\in N\left(i\right)} s_{ik}}}$$

(12)

In this study, MATLAB (MathWorks, USA) functions were used to implement the MFD algorithm. The function processes large DEM datasets and distributes water flow according to the MFD model, calculating the accumulated flow accurately for each pixel.

3. Results and Analysis

3.1. Centimeter-Scale Database of the Small Watershed

Based on UAV photogrammetry, a high-resolution (0.1 m) geographic information database of the small watershed was constructed (Figure 3). The elevation of the study area ranges from 56 to 119 m, exhibiting a typical north–south gradient characterized by higher elevations in the north and lower elevations in the south. The topographic gradient is pronounced across the watershed. The northern high-elevation areas are dominated by natural vegetation types, including bamboo forest, nursery land, grassland, and forest land, with high vegetation coverage and relatively stable ecosystem conditions. In contrast, the southern low-elevation areas are mainly occupied by intensively managed agricultural lands, such as vegetable fields, dryland, and greenhouse facilities, where human disturbance is strong.

A total of 14 land-use types were identified within the watershed. Natural land-use types constitute a substantial proportion of the area, with water bodies (26.09%) and forest land (15.26%) forming the ecological foundation of the watershed. Agricultural land-use types exhibit a complex spatial configuration. Dryland (10.37%), vegetable fields (5.01%), and greenhouses (2.96%) are clustered predominantly in the central part of the watershed, forming distinct zones of intensive agricultural activity.

Marked spatial contrasts were observed in both soil and streamwater nitrogen, indicating strong land-use controls and cumulative downstream inputs (Figure 4). Soil nitrogen varied substantially among land-use types. Total nitrogen (TN) was highest in greenhouse-based facility agriculture (2299.36 mg/kg), followed by vegetable fields and drylands (2006.68 and 1956.79 mg/kg, respectively), and was lower in paddy fields and bamboo forests (1742.36 and 1643.03 mg/kg, respectively). Ammonium nitrogen (NH₄⁺-N) peaked in drylands (39.86 mg/kg) and remained relatively high in greenhouse areas (37.98 mg/kg); vegetable fields, grasslands, and shrublands showed comparable concentrations (32.02–33.06 mg/kg), whereas paddy fields and bamboo forests had the lowest values (22.49 and 24.09 mg/kg, respectively). Nitrate nitrogen (NO₃⁻-N) was highest in greenhouse areas (73.65 mg/kg), relatively elevated in vegetable fields and grasslands (58.92 and 42.15 mg/kg, respectively), and lower in bamboo forests and paddy fields (48.97 and 40.88 mg/kg, respectively).

Streamwater nitrogen was dominated by downstream accumulation and responses to agricultural activities. Total nitrogen increased from 0.88–1.73 mg/L in the upper reaches to 1.38–2.44 mg/L downstream, frequently exceeding 2.0 mg/L near croplands and below greenhouse-dominated areas. NH₄⁺-N showed a similar pattern, remaining at 0.03–0.20 mg/L upstream but rising to 0.30–0.64 mg/L downstream; the agriculturally intensive area in the lower-right part of the watershed reached 0.58 mg/L. NO₃⁻-N increased from 0.22 mg/L upstream to a maximum of 0.64 mg/L downstream and was consistently higher in reaches downstream of cropland and greenhouse planting zones than in other areas.

3.2. Runoff and Sediment Transport in the Watershed

3.2.1. Runoff Simulation Based on the SCS-CN Model

The runoff process of the small watershed was quantitatively simulated using the SCS-CN model, and the results are shown in Figure 5. The spatial distribution of Curve Number (CN) values exhibits pronounced heterogeneity across the study area. Annual runoff shows clear spatial differentiation, with substantial differences in runoff generation capacity among land-use types. Paddy fields exhibit the highest runoff production, with an annual runoff depth of 530.1 mm. Roads, bare land, and residential areas also show relatively high runoff coefficients, resulting in runoff volumes significantly greater than those of natural vegetation areas. In contrast, agricultural land-use types such as vegetable fields, greenhouses, and dryland show relatively small differences in runoff depth. Bamboo forest, which has the highest vegetation coverage in the watershed, produces the lowest annual runoff among all land-use types. This spatial pattern of runoff provides a critical hydrological basis for subsequent simulations of pollutant loss.

3.2.2. Soil Erosion Estimation Based on the MUSLE Model

The spatial distributions of the MUSLE model factors are shown in Figure 6. High values of the soil erodibility factor (K) are mainly concentrated in cultivated land and bare land, whereas water bodies and forest land exhibit the lowest K values. Water bodies inherently have no erosion risk, while forest land shows low erodibility. The conservation practice factor (P) displays a spatial pattern positively correlated with the intensity of human activities, with higher P values observed in farmland and residential areas, indicating relatively poor soil conservation conditions; in contrast, bare land shows relatively lower P values. The vegetation cover and management factor (C) exhibits the most pronounced spatial variation, with values reaching 0.15–0.25 in bare land and greenhouse areas, while natural vegetation areas such as forest land and grassland show much lower values (0.01–0.06). The topographic factor (LS) is significantly higher in steep slope areas than in the relatively flat agricultural lands in the southwestern part of the watershed.

By integrating all MUSLE factors, the mean annual soil erosion intensity of the watershed is estimated to be 14.55 t/km². According to the Soil Erosion Classification and Grading Standard (SL190—2007), the watershed as a whole falls within the category of slight erosion. However, the spatial distribution of erosion intensity is highly uneven. Greenhouse areas, residential land, and bare land form distinct erosion hotspots. Among them, greenhouse areas exhibit the highest average erosion intensity, reaching 21.6 t/km², and represent the most severely eroded zones within the watershed.

3.3. Quantification of Nitrogen Output Loads in the Small Watershed

Based on the coupled simulation of the SCS-CN and MUSLE models, nitrogen output loads in different forms exhibit significant spatial heterogeneity across the small watershed (Figure 7). At the watershed scale, the mean annual output intensities are 29.66 kg/km²·a for total nitrogen (TN), 0.51 kg/km²·a for ammonium nitrogen (NH₄⁺-N), and 0.83 kg/km²·a for nitrate nitrogen (NO₃⁻-N). The spatial patterns indicate that hotspots of nitrogen export largely coincide with areas of elevated nitrogen concentrations and severe soil erosion, and are mainly clustered in the central greenhouse cultivation zone, drylands, and residential areas where human activities are intensive. In contrast, nitrogen export loads were markedly lower in areas dominated by natural vegetation cover, such as forests and grasslands. Notably, as shown in the figure, nitrogen export intensities along field boundaries were generally higher than those within field interiors, making these boundary zones critical nodes for non-point source pollutant export.

The regulatory effect of land-use type on nitrogen output loads is statistically significant, as confirmed by analysis of variance (ANOVA) (F = 8.96, p < 0.05). Quantitative results are summarized in Table 4. Bare land exhibits the highest total nitrogen output intensity (448.54 kg/km²·a), with ammonium nitrogen (10.51 kg/km²·a) contributing more than nitrate nitrogen (7.71 kg/km²·a). Greenhouse areas show a total nitrogen output of 363.41 kg/km², where nitrate nitrogen (12.89 kg/km²·a) exceeds ammonium nitrogen (6.41 kg/km²·a). Residential areas exhibit a total nitrogen output intensity of 239.29 kg/km². In contrast, natural ecosystems demonstrate strong nitrogen retention capacity. Total nitrogen output from forest land is only 6.09 kg/km²·a, with ammonium nitrogen at 0.13 kg/km²·a, while bamboo forest shows an extremely low total nitrogen output of 0.10 kg/km²·a, with negligible losses of both nitrate and ammonium nitrogen. Although paddy fields generate the highest runoff volume, their total nitrogen output intensity (6.68 kg/km²·a) is comparable to that of forest land.

Table 3. Mean output intensity of different nitrogen forms under different land-use types.

Land-Use Type	TN (kg/km²·a)	NO₃⁻-N (kg/km²·a)	NH₄⁺-N (kg/km²·a)
Bare land	448.54	7.71	10.51
Greenhouse	363.41	12.89	6.41
Residential area	239.29	5.64	3.35
Dryland	116.81	3.38	2.39
Vegetable field	93.35	2.66	1.49
Nursery land	40.74	1.02	0.73
Grassland	9.53	0.26	0.17
Shrub land	7.18	0.18	0.13
Paddy field	6.68	0.16	0.09
Forest land	6.09	0.17	0.13
Bamboo forest	0.10	0.00	0.00

Table 3. Mean output intensity of different nitrogen forms under different land-use types.

Land-Use Type	TN (kg/km²·a)	NO₃⁻-N (kg/km²·a)	NH₄⁺-N (kg/km²·a)
Bare land	448.54	7.71	10.51
Greenhouse	363.41	12.89	6.41
Residential area	239.29	5.64	3.35
Dryland	116.81	3.38	2.39
Vegetable field	93.35	2.66	1.49
Nursery land	40.74	1.02	0.73
Grassland	9.53	0.26	0.17
Shrub land	7.18	0.18	0.13
Paddy field	6.68	0.16	0.09
Forest land	6.09	0.17	0.13
Bamboo forest	0.10	0.00	0.00

3.4. Simulation of Nitrogen Migration Trajectories in the Small Watershed

Based on high-resolution topographic data and the multi-flow direction (MFD) algorithm, nitrogen migration pathways in the small watershed were simulated and visualized (Figure 8). Analysis of representative scenarios reveals the regulatory effects of microtopographic elements on nitrogen migration pathways. In scenarios where field ridges act as barriers between adjacent plots (Figure 8b), nitrogen migration exhibits a dual-path pattern: the majority of nitrogen is transported along ridge margins toward lower elevations, while a smaller proportion penetrates the ridges and enters neighboring fields. When ditches are present within fields (Figure 8d), nitrogen preferentially converges into the ditch network and is subsequently transported rapidly downstream through artificial drainage systems. In fields adjacent to water bodies (Figure 8c), nitrogen migrates directionally from cropland toward surface waters; however, part of the nitrogen load is intercepted within vegetated buffer zones. In areas surrounding buildings (Figure 8e), nitrogen migration pathways are significantly deflected by structural barriers, forming directional transport corridors along building edges.

The watershed-scale simulation results based on the MFD algorithm (Figure 9b) demonstrate pronounced spatial heterogeneity in nitrogen migration. Overall, nitrogen transport is closely related to terrain relief and the elevation gradient, forming a directional migration tendency from higher to lower areas. In the northern steep-slope areas, nitrogen migration exhibits typical linear convergence pathways along the slope direction, whereas in the southern low-gradient areas, migration transitions to a diffuse, planar spreading mode, with more dispersed transport pathways and a relatively slower migration process. Spatial differentiation in nitrogen accumulation is strongly consistent with land-use patterns (Figure 9a). In general, areas with higher vegetation cover, such as forest land, show lower accumulation levels, while agricultural land and residential areas with more intensive human activities are more likely to form accumulation hotspots.

Morphological differentiation analysis indicates that the overall migration extent is spatially constrained, with migration mainly concentrated in source areas and their immediate surroundings, and no cross-watershed-scale continuous convergence network is formed.

The source–sink pattern of nitrogen migration is strongly coupled with land-use structure (Figure 9g–j). Residential areas in the central-northern and northeastern parts of the watershed form the core source zones, while the convergence area at the junction of agricultural drainage ditches and surface waters in the southern part acts as the primary sink. This sink area collects multiple migration pathways from the northern residential zones and the central agricultural area, and serves as a key node for nitrogen spatial redistribution. Notably, greenhouse areas in the central watershed show a distinct dual role: on the one hand, they act as sink zones, receiving and accumulating nitrogen transported from upstream residential areas; on the other hand, they function as secondary source zones, continuously delivering both locally derived and accumulated nitrogen to downstream ponds and reservoirs through the ditch network.

4. Discussion

4.1. Interpretation of Field-Scale Nitrogen Migration Characteristics

The field scale represents a critical linkage between microscale soil processes and macroscale watershed responses. At this scale, nitrogen migration mechanisms are jointly regulated by topographic features and land-use types. Based on a 0.1 m high-resolution DEM, multi-model coupled simulations, and field observations, this study systematically elucidates the interactive mechanisms between topography and land use that govern nitrogen migration at the field scale.

Land-use type is the primary anthropogenic driver determining the intensity of nitrogen loss. Nitrogen outputs vary significantly among different land-use types. Specifically, greenhouse areas and vegetable fields exhibit markedly higher median nitrogen outputs than other land-use categories. This pattern can be attributed to intensive agricultural practices characterized by frequent and high-rate fertilizer applications and repeated soil tillage, which result in soil nitrogen surpluses that are readily mobilized by surface runoff [34]. In contrast, forest land and bamboo forest show the lowest nitrogen outputs, indicating that natural vegetation systems possess strong nitrogen retention and assimilation capacities. Although paddy fields involve artificial drainage processes, their nitrogen output levels are intermediate, which may be associated with denitrification processes in paddy wetland systems and optimized water–fertilizer management practices [33].

Topographic factors exert a fundamental influence on nitrogen migration through the combined effects of natural terrain properties and anthropogenic landscape modifications, shaping migration directionality, hydrodynamic behavior, and spatial redistribution processes. Consequently, topography serves as a core driving force controlling pollutant transport within agricultural landscapes [35].From the perspective of natural terrain, slope, slope length, and microtopography exhibit hierarchical control effects on nitrogen migration [36]. As slope increases, gravitational potential energy of surface runoff is converted into kinetic energy, significantly enhancing both the detachment and transport of soil particles and the mobilization efficiency of dissolved nitrogen [37]. This effect is particularly pronounced in agricultural environments such as dryland fields with low vegetation coverage, where slope-driven runoff substantially amplifies nitrogen loss risk. Slope length also exerts a positive influence on nitrogen output (Figure 11b), as longer slopes promote runoff accumulation by extending flow convergence pathways, thereby increasing nitrogen transport fluxes. However, comparison of coefficients of determination indicates that, at the field scale, slope explains nitrogen migration more effectively than slope length. This finding suggests that vertical gravitational forcing plays a more critical role in triggering nitrogen loss than horizontal flow convergence pathways.

From the perspective of anthropogenic microtopography, field ridges exhibit a decisive interception and regulation function. Some studies demonstrated that field ridges significantly disrupt lateral hydrological connectivity between adjacent fields, physically intercepting up to 91.1% of farmland runoff, with only approximately 10% of runoff generated at higher slope positions able to cross ridges[38]. The microtopographic simulations conducted in this study (Figure 8b) are consistent with these findings, further confirming the pivotal role of field ridges as artificial topographic barriers in modifying runoff pathways and mitigating cross-field nitrogen transfer.

4.2. Effects of DEM Resolution on Model Outputs

Nitrogen migration and transformation processes are tightly coupled with hydrological cycles and soil erosion dynamics. As the foundational topographic input for hydrological and soil erosion models, DEM resolution directly determines the accuracy of slope, slope length, and hydrological connectivity calculations. In agricultural small watersheds, microtopographic features such as field ridges and ditches constitute key structural elements controlling nitrogen source–sink transformation. The smoothing effect introduced by low-resolution DEMs can obscure these features, leading to misrepresentation of runoff pathways and sediment transport ratios and, consequently, substantial biases in nitrogen loss estimates. Therefore, quantifying the influence of DEM resolution on model performance is a prerequisite for ensuring the reliability of nitrogen migration simulations and for elucidating the regulatory role of microtopography in non-point source pollution control.

In this study, the original 0.1 m DEM was systematically resampled to seven resolution levels (0.5 m, 1 m, 2 m, 5 m, 10 m, 15 m, and 30 m) using ArcGIS 10.8, and soil erosion intensity under different resolutions was quantified using the MUSLE model (Figure 10). The results indicate that DEMs with resolutions of 0.1–0.3 m can accurately identify microtopographic features such as field ridges, furrows, and artificial drainage channels, with clear spatial structures. As resolution decreases below 5 m, these critical features become increasingly blurred, and at a resolution of 30 m, field boundaries and microtopographic details are almost completely lost. These findings corroborate the conclusion that low-resolution DEMs are inadequate for characterizing land-use-related microtopography[39].

Figure 10. Spatial distribution of soil erosion intensity at 0.1 m-30 m resolution DEM (a-h).

Figure 10. Spatial distribution of soil erosion intensity at 0.1 m-30 m resolution DEM (a-h).

Further analysis reveals a non-monotonic response of erosion intensity to DEM resolution (Figure S1a). The mean erosion intensity simulated at resolutions of 0.1–0.5 m (35.35 kg/km²·a) is significantly lower than that obtained at 1–2 m resolution (65.63 kg/km²·a). This phenomenon arises from the refined representation of microtopography in sub-meter DEMs, where local slopes are fragmented and micro-depressions and micro-elevations partially offset one another, resulting in an overall smoothing of slope gradients. As resolution continues to decrease, erosion intensity exhibits a pronounced inflection point at approximately 1 m resolution, followed by a declining and stabilizing trend [40].

To further elucidate this phenomenon, elevation differences of representative microtopographic features (e.g., field ridges and ditches) were measured in the study area. Field measurements indicate that elevation differences of typical microtopographic features are mainly distributed within the range of 0.5–2.0 m, with a median value of 1.2 m (Figure S1b). This characteristic explains the inflection point observed at 1 m resolution. At this resolution, terrain noise smaller than 1 m is effectively filtered out, while key topographic structures larger than 1 m are preserved, achieving an optimal balance between terrain detail representation and model stability [27]. Consequently, although sub-meter DEMs provide the most detailed terrain information, a resolution of approximately 1 m may offer greater practical advantages for field-scale studies in small watersheds.

Quantification of factor contributions based on structural equation modeling (SEM) further reveals the scale-dependent dominance of topographic controls on field-scale nitrogen migration (Figure 11). Under high-resolution scenarios (0.1 m, 0.5 m, 1 m, and 2 m), the contribution of topographic factors remains stable at around 80%, significantly exceeding that of vegetation cover (10–25%). This finding reinforces the critical role of microtopography in constraining nitrogen spatial redistribution processes. Topography influences nitrogen migration capacity by regulating flow velocity and kinetic energy through slope variation, while simultaneously determining migration pathways and accumulation patterns by controlling runoff convergence routes [41]. Together, these mechanisms form a dual topographic control on nitrogen migration.

When DEM resolution decreases to 5 m, the contribution of topographic factors declines to 40–50%, and at 30 m resolution, it further drops to 30%, indicating a loss of topographic control over migration pathways. In such cases, meteorological factors emerges as the dominant controlling factor. This shift is fundamentally driven by spatial attenuation of terrain information: low-resolution DEMs smooth micro-slope variations (e.g., elimination of slopes smaller than 2°), distortion of flow accumulation and direction allocation, and oversimplify hydrological connectivity networks, thereby weakening topographic regulation of nitrogen migration. Importantly, the increased contribution of meteorological factors does not reflect an actual enhancement of their controlling effect, but rather a compensatory dependence of the model on remaining explanatory variables following terrain information loss, resulting in a "pseudo-dominance" phenomenon [42].

At the field scale, the results of this study demonstrate that DEMs with resolutions of 0.1 m, 0.5 m, 1 m, and 2 m effectively preserve the spatial variability of small-scale terrain features such as ridges, furrows, and micro-depressions, thereby enabling accurate representation of topographic controls on slope-scale nitrogen migration. These microtopographic features play a crucial role in runoff pathway bifurcation and temporary water storage formation. High-resolution DEMs are capable of capturing these details and providing key predictive indicators for nitrogen migration, including terrain roughness indices and flow-path density. In contrast, low-resolution DEMs (30 m) fail to accurately delineate geomorphic unit boundaries, resulting in insufficient spatial precision and systematic bias in nitrogen load estimation [43].

Figure 11. Contribution of different factors under varying DEM resolutions.

Figure 11. Contribution of different factors under varying DEM resolutions.

5. Conclusions

This study, leveraging UAV photogrammetry and the coupled SCS-CN and MUSLE models, conducted a quantitative simulation of nitrogen migration characteristics and output loads at the field scale in a small watershed. The results are as follows:

High spatial heterogeneity of nitrogen output in the small watershed: The simulation results show that the mean annual total nitrogen (TN) output intensity in the watershed is 29.66 kg/km²·a. The highest output is observed in barren land (448.54 kg/km²·a) and greenhouse areas (363.41 kg/km²·a). High nitrogen load zones are concentrated in areas with intensive human activities and strong soil erosion, and field boundary areas, influenced by microtopography, serve as critical nodes for non-point source pollution output.

Nitrogen migration trajectories are predominantly influenced by topographic gradients, with significant influence from the synergy between topography and land-use types at the microscale: At the macro scale, nitrogen migration in the northern steep-slope areas exhibits a linear convergence pattern, while the southern gentle-slope areas exhibit a diffuse, planar spreading mode. At the microscale, land-use types determine the potential intensity of pollution sources, while microtopographic features, such as field ridges, ditches, and buildings, physically block or direct nitrogen migration, thereby altering the spatial distribution of nitrogen within different land-use units.

The resolution of DEM significantly affects the accuracy of nitrogen migration modeling and the identification of key controlling factors: The simulation results indicate a critical scale effect near 1 m resolution. Erosion intensity shows a non-monotonic relationship with resolution, with a distinct inflection point at 1 m, corresponding to the median elevation difference (1.2 m) of typical microtopographic features (e.g., field ridges) in the study area. This suggests that, under the geomorphic conditions of this study, the effective resolution for identifying key topographic features lies around 1 m. Structural equation modeling further confirms that high-resolution (0.1–2 m) DEMs effectively capture microtopographic features, maintaining the contribution of topographic factors to nitrogen migration around 80%. However, when the resolution drops to 30 m, the contribution of topography significantly decreases to 30%, resulting in a marked shift in the identification of dominant controlling factors.