Influenced by solar radiation and gravity, various dissolved and insoluble substances (such as salts, gas molecules and sediments) migrate continuously between the atmosphere-soil, atmosphere-water and soil-water interfaces through several hydrological paths, such as precipitation, infiltration and runoff (Taikan et al., 2006; Yang et al., 2011; Pang et al., 2012; Masese et al., 2016; Ning et al., 2019). During these processes, the amount of mass that flows through such paths per unit area per unit time under conditions of hydraulic circulation is the biogeochemical mass interface flux (hereinafter referred to simply as “flux”) (Kong et al., 2015; Miao et al., 2016; Su et al., 2019). The magnitude of flux can reflect a considerable amount of information that is closely associated with biogeochemical basin processes (Del Vecchio et al., 2018; Moore et al., 2008; Nishino et al., 2018; Lambrecht et al., 2019), which could also even directly affect primary productivity (Figure 1). For example, increased atmospheric nitrogen (N) deposition may result in N saturation within a watershed ecosystem (Gao et al., 2019, 2020), which is defined as the condition wherein N availability exceeds the capacity of plant and soil microbial requirements, which can lead to noxious and toxic algal blooms, increased turbidity, disruptions to ecosystem functions, shifts in the food web and losses in fish stocks (Gundersen et al., 1998; Rabalais, 2002; Dentener et al., 2006). Moreover, acid deposition can increase dissolved organic carbon (DOC) in surface water of glaciated landscapes, which can subsequently increase the export of DOC to ocean systems, a potentially important component of the regional carbon (C) balance (Jan, 2003; Monteith et al., 2007). Land-use change will affect mass flux in the mid-layers of soil, which is an important means of soil nutrient loss, significantly impacting plant nutrient absorption (Bhaduri et al., 2000; Turner et al., 2001; Zhu et al., 2014). Furthermore, sediments process a considerable amount of nutrients that are input into systems, which facilitates both ecological productivity and ecological diversity (Battin et al., 2016; Voermans et al., 2018). In addition, the fluxes of other micromolecules (such as phosphate, metal ions, etc.) and macromolecules (such as black carbon, etc.) are also important components of the biogeochemical cycle in watersheds (Coppola et al., 2018; Gao et al., 2018).

In practical applications, the real value of flux is unknown; thus, it can only be estimated. The accuracy and precision of flux estimation results are affected by measurement errors and system errors. Measurement errors refer to errors caused by monitoring methods, such as velocity measurement errors, water quality sampling errors, water quality analysis errors, discrete cross-sectional sampling errors, sampling frequency errors, etc. (Zhao et al., 2017; Zhao et al., 2018). Systematic errors refer to errors caused by the utilization of different estimation methods. The accuracy and precision of flux estimation results directly affect the assessment of environmental effects and the perception of the environmental status, ultimately affecting the formulation and implementation of watershed management planning. Therefore, accurately estimating flux is important in watershed ecology (Baldassarre and Montanari 2009) and provides a scientific basis for the sustainable development of both regional and global environments. Scientifically grasping the law of variation associated with pollution loads during mass-hydrological coupling processes can provide a scientific basis for the study of environmental capacities as well as the control and reduction of total pollution loads (Lovett, 1994; Zhao et al., 2014).

It is important to understand that flux obtained from different monitoring and estimation methods will significantly differ (Webb et al., 1997). However, existing studies on mass flux associated with biogeochemical watershed interfaces lack uniform criteria or references as it pertains to the application of monitoring and estimation methods. In order to resolve this issue, this study summarized and deconstructed results from recent studies on biogeochemical watershed interface processes and compared the advantages, disadvantages and applicability of the monitoring and estimation methods used by various relevant studies. This particular study is intended to be used as a reference for the selection of flux calculation methods.

Watershed water cycling processes include deposition, soil leaching, river export and evapotranspiration, etc. Numerous dissolved or insoluble substances can migrate and transform under these processes (except for evapotranspiration). The sediment-water interface (SWI) is an important interface for mass exchange between solids and liquids on the Earth’s surface. Therefore, this study focused on four different types of mass transport flux processes, namely, deposition, soil leaching, river export and SWI diffusion (Figure 1).

In theory, the following formula is typically used to estimate the flux of a substance over a period of time:

where W is flux (m³); Q(t) is the instantaneous flow rate (m³/s); C(T) is the instantaneous concentration (mg/L) (Wang et al., 2004; Armstrong et al., 2010); C_i is the instantaneous concentration (mg/L); Q_i is the instantaneous flow rate (m³/s); ∆t is the monitoring interval time; n is the monitoring time (Allan et al., 2006; Hao et al., 2017; Lu et al., 2019). Flux values obtained using different estimation methods differ significantly; thus, the estimation method used should be carefully considered (Liu et al., 2019). In general, flux estimation should be as precise as possible on the basis of accuracy (Rekolainen et al., 1991; Zhang et al., 2015).

Deposition is the main transportation pathway of airborne contaminants, such as trace metals, nitrates, sulfates and toxic pollutants, moving from the atmosphere to terrestrial and aquatic ecosystems (Kim et al., 2001; Landis and Keeler, 2002; Rolfhus et al., 2003; Huang et al., 2011). Atmospheric deposition includes dry deposition and wet deposition (Ottley and Harrison, 1991; Andersen and Hovmand, 1999; Hayashi et al., 2013) (Figure 1). Dry deposition refers to the process wherein particles in the atmosphere are adsorbed into the terrestrial surface under forces such as gravity when precipitation does not occur. Wet deposition refers to the process wherein water droplets or ice crystals in the upper atmosphere carry airborne particles to the terrestrial surface when precipitation events occur (rain, snow, fog, etc.). At present, an international wet deposition observation network has basically been established. The methods used for wet deposition are a precipitation collection method and an ion exchange resin method, both of which have had time to mature. Moreover, dry deposition is more difficult to measure than wet deposition; thus, the construction of a dry deposition observation network is still in its infancy (Huang et al., 2011). Commonly used methods for obtaining dry deposition flux are the dust collecting cylinder method and the model simulation method.

2.1.1 Monitoring and estimating wet deposition flux

(1) Precipitation collection method

The precipitation collection method refers to the method of collecting precipitation regularly and measuring mass concentrations in rainwater. It can be divided into the “manual” collection and the “automatic” collection of precipitation. During the manual collection process, stainless steel instruments or polyethylene plastic barrels are typically used to collect samples. If rainwater is collected continuously over a long period of time, it is necessary to add antifungal agents or to rapidly secure samples in a refrigerator to prevent the transformation of different N forms during sample collection. The advantages of the artificial precipitation collection method are its low cost and the flexible sample arrangement it offers. The disadvantages of the artificial precipitation collection method are its higher observational requirements and its time-consuming and laborious nature. Therefore, this method is more suitable for use within field ecological stations under strong technical support (Sheng et al., 2010). When calculating deposition flux by the precipitation collection method, the following formula is used (Holland et al., 2005):

where R_w is the wet deposition rate (kg/km/month); C_i is the concentration of the i-th rain component (mg/L); h is rainfall (mm); k is the unit conversion factor.

The automatic precipitation and dust collector methods employing automatic precipitation sample collectors have emerged in recent years. When a precipitation event occurs, the instrument, under the control of a rain sensor, will automatically collect rainwater while closing the dust collection apparatus (Holland et al., 2005; Xi and Mulder, 2007). The advantages of these precipitation and dust collectors are that they can be operated by a microcomputer, and thus record precipitation automatically; certain new types of precipitation collectors are equipped with refrigeration equipment to ensure that precipitation samples do not deteriorate during sampling. The disadvantage of these methods is that the precipitation and dust collectors require a stable power supply; thus, these collectors are difficult to maintain and apply in remote locations.

(2) Ion exchange resin method

Functional groups in ion exchange resins can dissociate certain cations (e.g., H⁺ or Na⁺) or anions (e.g., OH^-or Cl^-) in aqueous solutions, while adsorbing other original cations or anions in solutions. Through ion exchange, the ions that are to be measured in precipitation are fixed on the functional groups of resins (Klopatek et al., 2006; Sheng et al., 2013). The advantages of this method are that requirements for sample preservation are low, which facilitates the elimination of human sampling errors between stations; and cloud deposition can be captured, which provides more realistic results from forest ecosystems in warm wetland zones compared to those obtained by traditional methods. The disadvantages of this method are that the higher sampling temperature requirement and the potential influence of adsorption and desorption processes will affect test result accuracy; and test result accuracy will be significantly influenced by progressively aging resin (Fenn and Poth, 2004). The formula to calculate deposition flux using the ion resin method is as follows (Sheng et al., 2010):

where D_IER is the sedimentation flux (kg/hm); C_ex is the content of the extract (mg/L); V_ex is the volume of the extract (L); A is the area of the funnel (m³); 100 is the unit conversion factor.

2.1.2 Monitoring and estimating dry deposition flux

(1) Wet collection method of the dust collector

The wet collecting method of the dust collector uses an organic glass dust collecting cylinder to collect dry sediment samples from the atmosphere. Dust collectors are typically placed at a height of 1.2 m relative to the surface away from large trees or buildings. Additionally, 5 cm of distilled water should be added to the dust collector. The cover must be sealed against precipitation, while dust samples must be collected following rain events. The advantages of using the wet collecting dustfall cylinder method are its simple operation and its low cost. It is therefore one of the most common methods used to monitor atmospheric dry deposition (Chen et al., 2006). The disadvantages of using the wet collection dustfall cylinder method are the inconvenience of adding water into the dust collector; the excessive rainfall during the wet season makes it easy for water to overflow into the dust collector; and water in the cylinder can easily evaporate in summer and easily freeze in winter. During field operations, this method is often substituted for automatic precipitation and dust collectors to allow for automated switching between wet and dry deposition collection processes.

(2) Model simulation method

The model simulation method is an alternative method that is mainly used to calculate atmospheric dry deposition flux in a large monitoring network (Hicks et al., 1986, 1987). For this method, atmospheric deposition flux is equal to the product of the concentration of mass and deposition rates (Baker et al., 2010). Among these, mass concentrations in the atmosphere are measured using special instruments, while the sedimentation rate is calculated by the model after measuring meteorological parameters and underlying surface parameters (Wang and Li, 1994; Mao and Hu, 1996; Pratt et al., 1996; Yan et al., 2004). Due to its simplicity, this method has been widely used to monitor and calculate dry deposition flux (e.g., CASTNet in the United States and CAPMN in Canada) (Fan et al., 2007, 2009). The formula used to calculate dry deposition flux of atmospheric particulates under model simulations is as follows (Wesely et al., 2000; Zhang et al., 2001, 2003; Holland et al., 2005; Lyman et al., 2007):

where F_c is the dry deposition flux of atmospheric particulates; C_a is the average concentration of particulates in gases and aerosol particles; V_d is the dry deposition rate. Deposition velocity is assumed to be inversely proportional to the sum of three resistance factors:

where r_a is the aerodynamic resistance; r_b is the resistance of the laminar sublayer between the surface and the turbulent boundary layer, which depends on the diffusivity of the species; r_c is the resistance of the bulk surface characteristics and their correlation to the solubility and reactivity of chemical species. The advantages of the simulation method are the low sensitivity requirement of sensors; the relatively simple measurement process, suitable for long-term monitoring; and that it can be combined with a geographic information system (GIS), after which estimations of regional dry deposition flux can be realized.

(3) Other methods

In addition, there are other albeit imperfect methods to estimate dry deposition flux, such as canopy budget models for forests based on throughfall and stemflow measurements (Hayashi et al., 2013; Aguillaume et al., 2017) and the passive dry deposition (Pas-DD) collector. The Pas-DD collector is promoted as a tool to estimate atmospheric loadings (dry deposition flux or maps) for polycyclic aromatic hydrocarbons (PAHs) and other semivolatile compounds (Eng et al., 2013).

Soil leaching refers to the process of vertical (to the depth of the soil) or horizontal migration of soluble substances or suspended compounds in the soil along with interflow. Interflow refers to the movement of water in soil, which includes vertical infiltration and lateral flow (Figure 1). During any rainfall event, at least a portion or potentially all water will infiltrate into soil layers along soil pores to form soil water. Scholars have conducted extensive research on the generation mechanism (Beven 1989), preferential path (Zhang et al., 2004; Uchida et al., 2005), criticality and nonlinearity (Tromp-Van Meerveld and Mcdonnell, 2006), and the influencing factors (Kienzler and Naef, 2008) of interflow.

2.2.1 Methods used for soil solution collection

Soil solution collection methods include both destructive and non-destructive collection.

(1) Destructive collection

Destructive sampling requires in-situ removal of soil samples, a process which includes centrifugation, the extraction method, the displacement column method and the pressure filtration method. Centrifugation is a method of separating soil solutions from soil applying high-speed centrifugal force. Its advantages are that soil solutions can correspond to specific soil layers in succession, and that soil water characteristics will not change. The disadvantages are that both the chemical composition and equilibrium of soil solutions are prone to change; it is difficult to conduct long-term positioning research; and it may not be possible to obtain sufficient soil solutions from soil with a low water content. The extraction method requires mixing and oscillating soil samples with water or diluted salt solutions at certain proportions. After water and soil are separated by filtration, centrifugation or dialysis, the ion concentration in the solution is determined. This method is suitable for investigating soil equilibrium processes, such as ion exchange, dissolution and precipitation. The replacement column method uses an effluent to wash fresh soil samples in the soil column, after which it is used to determine the ion concentration in the effluent. The ion composition of the soil solution could however change via the effluent reaction of certain ions. Pressure filtration is the process of removing the soil solution from soil by replacing the effluent with air pressure (Lawrence and David, 1996). The advantage of the pressure filter method is that sampling can achieve centimeter-level accuracy, while the disadvantages of this method are the maintenance of the pressure filter and that cleaning is time-consuming, requiring approximately 1 hour per sample (Böttcher et al., 1997).

(2) Non-destructive collection

The non-destructive collection method can be used for long-term soil positioning, such as measuring dynamic changes in soil solutions during plant growth, including the osmometer method, the negative pressure method, the diffusion method and the capillary method. For the osmometer method, the soil solution concentration is measured by a non-pressure osmometer, which collects soil moisture (i.e., gravity water) migrating downward along the soil substrate under the force of gravity (Giesler et al., 2010). This method is suitable for investigating the balance between inputs and outputs in ecosystems (Marques et al., 2010). However, the installation of non-pressure osmometers greatly disturbs soil (Giesler et al., 2010). The widely used negative pressure method (i.e., the suction cup method) is a standard method used to characterize hazardous waste points stipulated by the United States Environmental Protection Agency (EPA) (Brandi-Dohrn et al., 1996). When capillary pressure of the suction cup is less than that of soil, soil water is sucked into the suction cup (Grossmann and Udluft, 2006). Unlike non-pressure osmometers, suction cups can collect both gravity water and a portion of capillary water (Marques et al., 2010). The capillary method is associated with micro-analysis technology developed in recent years with the emergence of the micro-soil solution sampler and the application of micro-analysis technology. Its detection range is comparable to that of ion chromatography, but sample consumption is lower (only 5-10 nL per detection) and separation speed is faster (Wiltshire et al., 1995; Brandi-Dohrn et al., 1996). The diffusion method is used to collect soil solutions applying the passive diffusion principle conceived by Moutonnet (Moutonnet et al., 1993). For the method, the ceramic cup of a sampler contains degassed deionized water. After a balance is achieved at approximately 6-10 days, the soil solution reaches equilibrium with the solution in the ceramic cup (Moutonnet et al., 1993; Moutonnet and Fardeau, 1997).

Coinciding with the application of new technologies and new substances, new soil solution sampling techniques have and will continue to emerge. For example, a composite probe can simultaneously measure soil tension and collect soil solutions (Baumgartner, Parkin, and Elrick, 1994). In addition, the application of a time domain reflectometer has made it possible to monitor soil moisture and solute distributions rapidly, continuously and automatically at multiple points (Wraith et al., 1998). However, soil solutions collected by these different methods cannot be compared together owing to the different sampling devices, sampling principles and hydrogeological conditions of the sampling points used (Marques et al., 2010). Therefore, when choosing a sampling method, research objectives, accuracy requirements and characteristics of the research objects must be taken into account.

2.2.2 Methods used to monitor and estimate soil leaching flux

Calculation methods used for soil leaching flux include the in-situ plot method, the soil tank simulation test method and the model simulation method.

(1) In-situ plot method

Currently, the in-situ plot method is the most intuitive and accurate method to study slope runoff (Xiao et al., 2016). The runoff plot experiment method can simultaneously monitor both the flow and water quality of rainfall runoff and estimate flux within a study area as a whole to flux within a plot. Operation methods are as follows: selecting a typical representative runoff plot within a study area, excavating a section underneath it, building a surface runoff catchment trough at the foot of a slope to collect surface runoff and sediment, building a trench along the rock-soil interface to collect interflow and, finally, collecting lateral surface seepage at the bottom of a slope section. A soil moisture sensor and a soil water potential sensor are used to monitor dynamic changes in water content and water potential in the soil profile, respectively (Gao et al., 2014). The advantages of the in-situ plot method are its small workload and its cost effectiveness. Its disadvantages are the difficulty in determining typical plots; low flux estimation accuracy; and the understanding of regional pollution differences by only substituting large areas with plot studies is not conducive to the comprehension of regional pollution differences (Chen et al., 2014):

where Q_i is the soil leaching flux (mg/m²); c_i is the concentration of a substance in interflow (mg/L); q_i is the flow depth of the interflow (mm).

(2) Soil tank simulation test method

The soil tank simulation test method is used to simulate and calculate soil leaching through laboratory tests (Lei et al., 1988; Gu, 2000; Wang et al., 2012). Operation methods are as follows: A rectangular poly (methyl methacrylate) (PMMA) seepage experimental tank with a length of 50 cm, a width of 5 cm and a height of 20 cm should first be constructed, wherein the tank body is divided into three sections: a water chamber section (10 cm long, connected to a Marriotte bottle), a filter section (10 cm long, filled with fine quartz) and a sample section (30 cm long filled with soil samples). The outlet of the other end of the experimental tank is connected to a 1 L fine mouth bottle to receive seepage. The advantages of the soil tank simulation test method are that it is easy to control experimental conditions and observe experimental results, and that there is no need for long-term field observations. The disadvantage of this method is that it is impossible to exclude potential deviations between laboratory test conditions and natural conditions. The formulas to calculate soil leaching flux using this method are as follows:

where L_a is the leaching flux per unit mass of soil (mg/kg); L_p is the concentration of the seepage fluid (mg/L); V is the volume of the seepage fluid (L) calculated from rainfall; W is the mass of the soil samples (kg).

where L_s is the soil leaching flux per unit area (kg/hm²/a); ρ_b is the soil density (g/cm³); V is the total volume of soil per tillage layer (2×10³ m³).

(3) Model simulation method

The estimation of soil leaching flux using the model simulation method requires combining measured soil solution concentrations (obtained by the method discussed in section 3.2.1) with model estimated soil water flow, through means such as MIKE SHE (Demetriou and Punthakey 1998), MODFLOW (Tonkin et al., 2003), CENTURY (Liu et al., 2000), SOILN (Johnsson et al., 1987), DNDC (Li et al., 2006), HYDRUS (Lai et al., 2016; Šimůnek et al., 2016) and SWAP (de Vries et al., 2010) as well as other river basin hydrological models. The mathematical and physical methods used to determine the occurrence and movement of interflow are mainly based on the capillary potential theory proposed by Buckingham and the unsaturated flow equation proposed by Richards (Xiao et al., 2016). According to the description of interactions between soil and hydrological processes, these models are classified into three categories: a water storage-discharge model, a dynamic wave model and the Richards mode (Sloan and Moore, 1984; Li and Pei, 1999). The model simulation method typically requires additional watershed parameters. For example, to estimate soil water flux, the SWAP model requires meteorological parameters, abiotic parameters (such as drainage characteristics, physical soil characteristics, etc.) and biological parameters (such as the leaf area index, tree height, root index, etc.).

The accurate estimation of nutrient loads in rivers and streams is critical for many applications (Stenback et al., 2011). The instantaneous flux of river substances is the product of instantaneous flux and concentrations; however, to obtain long-term flux (such as annual flux), it is necessary to integrate monitoring records of instantaneous flux and concentrations for a specified period, while river basin runoff devoid of observational data can only be calculated by modeling (Cao, 2015; Gao et al., 2015). The accuracy and precision of river output flux estimations are influenced by the sampling frequency, sampling method, calculation method, watershed size and the type of chemical substances monitored (Richards and Holloway, 1987; Johnes, 2007; Hao et al., 2012). Compared to large watersheds, output processes of substances within small watersheds are closely related to underlying surface conditions, rainfall processes and anthropogenic activities. Problems associated with this procedure are considerable lag, fuzziness, latency, randomness and time differences (Li and Li, 2008; Kovács et al., 2012; Duan et al., 2014).

Determining the selection of the smallest sampling frequency to obtain the most representative data while minimizing flux estimation errors has always been a hot topic in the study of riverine material flux. The cost of collecting and analyzing water quality samples is high; thus, water quality data are typically collected at a lower frequency (Brauer et al., 2009; Duan et al., 2014; Park and Engel, 2014; Park et al., 2015), usually from 8 to 24 times per year for periodic and discontinuous monitoring (Allan et al., 2006; Ullrich and Volk, 2010; Hao et al., 2012). Some studies select sampling frequencies randomly and often provide no statistical basis for their results (Strobl and Robillard, 2008). Other studies search for a suitable monitoring sampling frequency to establish a widely applicable water quality monitoring program. However, due to the inherent differences in research conditions, different results are often obtained. For example, Kronvang and Bruhn (1996) studied two small watersheds in Denmark and determined that a sampling frequency of once every two weeks was adequate to monitor total N, total phosphorus (P), particulate P and soluble P (Kronvang and Bruhn, 1996). On the other hand, Coynel (2004) reported that within an error range of 20%, the monitoring frequency of suspended particulate matter should at least be once every 3 days for large watersheds and at least once every 7 hours for small watersheds (Coynel et al., 2004). To monitor N flux in small agricultural watersheds, local conditions and monitoring costs must also be taken into account (Valiela and Whitfield, 1989). Sample volume can be appropriately increased at the peak of flow to improve the accuracy and precision of sampling (Rekolainen et al., 1991).

Calculation methods used for river output flux mainly include the concentration-flow method, the empirical model method (based on empirical equations) and the mechanistic model method (based on physical mechanisms). Advantages of the concentration-flow method are ease of use and ease of calculation. Its disadvantages are high spatial and temporal data requirements and the high cost of long-term large-scale monitoring. The advantage of the empirical model method is its low data volume demand, while its disadvantage is its less accurate simulation results. Mechanistic models are mostly used to simulate hydrological processes or the migration and transformation of nutrients. The advantage of the mechanistic model method is that it can provide more accurate simulation results. Its disadvantage is that a large number of measured input parameters are needed to calibrate the model; thus, it is to a certain extent limited by the availability of parameters (Singh et al., 2005).

2.3.1 Concentration-flow method

(1) Method to monitor concentrations

Due to the uneven distribution of river water concentrations in monitoring sections, the adoption of a method which averages multiple sampling points is advantageous. For example, three sampling lines (left, middle and right) are arranged for each sampling section, and three sampling points are selected for each vertical line. The points on each vertical line are respectively 0.5 m below the water surface, the midpoint of water depth and 0.5 m above the river bed (Zhang et al., 2013). The substance flux monitoring section should be placed adjacent to the estuary (to control the entire basin area) while not being positioned by recharged water (Liu et al., 2011; Zhang et al., 2013).

(2) Method for flow monitoring

Flow monitoring methods include the buoy method, the flowmeter method, the volume method and the overflow weir method (Zhu et al., 2008). The buoy method is a simple method to roughly measure flow velocity. When measuring, a typical section of a river bed of a length no less than 10 m and of a specified liquid level height is selected. A floating substance is taken and placed in the middle of the upper reaches of the river. Under conditions where there is no external force (such as wind, blockage from floating substances, etc.), flow time of the float logistics is recorded via a measured distance. This is repeated several times to ascertain the average value. Flow is then calculated according to the following formula:

where Q₁ is river flow calculated by the buoy method (m³/s); L is the selected channel length (m); t is the average time (s) required for the buoy method; S is river section area (m²). The flowmeter method is suitable for measuring rivers with a water depth greater than 0.05 m and a flow velocity greater than 0.015 m/s. When measuring, the river section is typically divided into several sections, and the area and average velocity of each section are calculated separately, while river flow is calculated using the following formula:

where Q₂ is river flow calculated by the flow meter method (m³/s); F_n is the cross-sectional area of a section (m²); $\;{\bar v_n}$is the average flow velocity (m/s) of river water within a section. The volume method is used to import river water into a container of a known volume, measure the time it takes to fill the container and use this data to divide the volume of the receiving container to obtain the flow rate. For the overflow weir method, triangular, rectangular or trapezoidal weir plates are used to block water flow to form overflow weirs (Figure 2 represents a right-angled triangular weir). Flow calculation formulas are as follows:

where Q₃ is river flow calculated by the volume method (m³/s); H is head height (m); K is the discharge coefficient; D is the height (m) from the bottom to the edge; B is the width (m) of the upper flow of the weir.

(3) Method to calculate flux

Calculation methods for flux estimation using the concentration-flow method are mainly subdivided into four types: the sum of time-interval flux, the product of time- interval average concentration and water volume, the sum of flux frequency distribution and the convection-diffusion model. Among these, the first two types are more commonly used, while the second type is less accurate than the first type and the last type is only applicable to branching estuaries. Webb (1997) constructed five time-interval flux estimation formulas based on the first two calculation types (Table 1). Methods A and B only contain products of mean flow rates and mean concentrations, neglecting discrete time-averaged terms, while methods C, D and E include both. Method F is used to establish logarithmic relationships between concentration and the flow rate, while method G applies the correction coefficient (CF2) on the basis of method F (Webb et al., 1997; Johnes, 2007). Although both methods A and B underestimate suspended sediment flux, they are relatively accurate. Method B overestimates dissolved chemical flux under a simple dilution response. The estimation of suspended sediment flux using method F is relatively low (Hao et al., 2012; Zhang et al., 2015). Method C is suitable for pollutants whose flux is not closely associated to the flow rate. Methods D and E (wherein flow rate correction causes deformation in method D) are both suitable for pollutants whose flux is closely associated with the flow rate (Hao et al., 2012). Method C is superior to the other estimation methods in most cases, and this is simply because most rivers are polluted by both point and non-point pollution sources (Wang et al., 2011).

2.3.2 Empirical model method(1) General empirical formulaWhen measured data are scarce, the flux of some specific substances can be estimated using empirical formulas. Researchers have proposed many empirical models under different application ranges. For example, the empirical formula for calculating DOC flux was proposed by Ludwing:

where F_DOC is DOC flux; Q is runoff depth (mm); S is the slope of the topography (unit: radian); C_soil is the organic soil C content (%) (Ludwig et al., 1996; Gao and Shen, 1998). For another example, Schlesinger’s empirical formula for calculating organic C flux is as follows:

where Y is organic C flux (g); X₁ is total annual runoff (m³); X₂ is watershed area (km²) (Schlesinger and Melack, 1981; Gao and Shen, 1998).(2) Box model methodThe box model can be used to estimate mass flux in river estuaries (Worrall and Burt, 1999; Webster et al., 2000; Sun et al., 2006). The principle of the box model method is to divide flow data into several boxes according to corresponding time intervals. Assuming that the flow rate and the concentration of the particle in each box are the same according to the characteristics of flow data corresponding to the sample, the value of the sample is distributed to all boxes.(3) Rating curveThe most common method to combine intermittent concentration data with continuous discharge data uses a rating curve to predict unmeasured concentrations from discharge at a specific time; however, most estimate results using this method will be too low (Ferguson, 1986).The rating curve almost always uses the least squares regression equation:

where logarithms are base 10, and a and b are constants.(4) Load duration curveBased on the flow duration curve (FDC), the load duration curve (LDC) is the relationship curve between mass flux and duration range. LDC can directly reflect substance flux change characteristics and accurately estimate annual river flux under large annual flow rate variation. LDC is easy to understand and is extremely suitable for riverine applications where hydrological and water quality monitoring data are relatively scarce. LDC uses two flux estimation methods (Horvath et al., 1998; Li et al., 2012):① The average concentration estimation method:

where F_c is the substance flux; $\bar Q$ is the average flow rate within a divided time interval; $\bar C$ is the average concentration of a particle (mg/L); N is the total number of time intervals; n is days within the time interval.② The average flux estimation method:

where F_f is the substance flux; Q is the daily flow rate; $\bar C$ is the average concentration of a particle (mg/L); $Q\bar C$ is the average daily flux.(5) Export coefficient modelThe Export Coefficient Model (ECM) is widely considered to be a reliable method for simulating non-point source pollution loads in large and medium-scale watersheds (Johnes, 1996; Worrall and Burt, 1999). This method makes full use of available land-use data and directly establishes the relationship between the non-point source pollution load and the land-use type (Geng et al., 2013). Due to the wide application of this model, many researchers have improved upon it and have consequently considered pollution sources more exhaustively. Among these improved ECMs, the one proposed by John (Johnes, 1996) is the most representative. This improved ECM (Ding et al., 2010) reduces the influence of rainfall and runoff on simulation results. This formula is as follows:

where L is the amount of mass loss (kg); E_i is the output coefficient of the i-th nutrient source (kg/km²/a); A_i is the area of the i-type land-use type (km²) or the number of the i-th livestock; I_i is the nutrient input amount of the i-type (kg); P is the mass of the substance input through precipitation (kg); α is the rainfall impact factor; β is the topographic impact factor.(6) LOADEST modelThe LOADEST model uses continuous water volume data and discrete water quality data in establishing the mass flux regression equation (Park and Engel, 2014). A web-based tool that uses LOADEST as a core engine supporting four modules was developed to provide user-friendly interfaces and input data collection via web access (Spencer et al., 2013; Pellerin et al., 2014; Park et al., 2015). In addition, the LOADEST model can also be used to determine the optimal frequency of water quality monitoring (Li et al., 2012; Zhang and Chen, 2014). Examples of similar regression-based models include the composite method (Aulenbach and Hooper, 2006) and WRTDS (Hirsch et al., 2010).(7) SCS modelThe USDA runoff curve number method (commonly referred to as the SCS model) has proven to be an enduring method for estimating the volume (and peak rate) of direct (surface) runoff from ungauged rural catchments (Liu et al., 2011; Xu et al., 2018). It came into common use under the auspice of the United States Department of Agriculture (USDA) Soil Conservation Service (SCS) around the mid-1950s and continues to be used as the basic hydrological model that is incorporated into many USDA SCS systems, such as CREAMS and SPUR (Boughton, 1989). The model used for estimating surface runoff was also developed by the USDA SCS. It has been widely used to estimate storm runoff depth in all patches within a watershed based on runoff curve numbers (CN) (USDA, 1972). The SCS equation for storm runoff depth is mathematically expressed as (Weng, 2014):

where Q is storm runoff; P is rainfall; S is the potential maximum storage capacity, and $S = \frac{{1000}}{{CN}}$; CN is the runoff curve number of the hydrologic soil group-land-cover complex.2.3.3 Mechanistic model methodThere are many types of mechanistic models, such as ANSWERS, SWRRB, SWAT, HSPF, AGNPS, BASINS, SWMM, STORM, SLAMM, etc. (Liao et al., 2013). Compared with the empirical model method, the mechanistic model method has higher accuracy, but requires more input parameters. Below are the four main types:(1) SWAT modelThe Soil & Water Assessment Tool (SWAT), being a widely used model in simulating non-point source pollution in river basins, can simulate and calculate runoff, soil erosion and nutrient loads over a long period of time. It has become an indispensable tool in water resource protection and management (Grizzetti et al., 2003; Wang et al., 2008; Chen et al., 2010; Ullrich and Volk, 2010a; Geng et al., 2013; Tobin and Bennett, 2013; Gui et al., 2014; Li et al., 2016).(2) HydroTrend modelBased on the water balance principle, taking into account both topographic and geological basin characteristics (i.e., elevation, river networks, slope and lithology), physical parameters (i.e., temperature, rainfall, evapotranspiration, ice and snow cover) and anthropogenic activities (i.e., soil reclamation and dam construction) HydroTrend 3.0.4 is a lumped hydrological model used to simulate water and sediment flux at river basin outlets (Syvitski et al., 1998; Kettner and Syvitski, 2008; Sheng et al., 2018).(3) The Mike Basin modelThe Mike Basin is an ArcGIS-based integrated planning and management tool used for water resources on a watershed or regional scale. This mathematically-based software can simulate surface water runoff (Cao, 2015) as well as river material output flux.(4) RHESSys model (Son et al., 2019)The Regional Hydro-Ecological Simulation System (RHESSys) is a biophysically based, spatially explicit model that has the capacity to simulate both climate and vegetation change impacts and feedbacks between hydrologic processes and ecosystem C and nutrient cycling.

The sediment-overlying water interface is an important interface for particle circulation in lakes and reservoirs (Pitkänen et al., 2001; Mu et al., 2017; Voermans et al., 2018; Zhao et al., 2018). Understanding the transport processes of the SWI is fundamental in understanding the role that sediments play in aquatic ecosystem functions (Huettel et al., 2014; Voermans et al., 2018). The exchange of particles at the SWI is a two-way process (i.e., where particles entering sediments through water and water through sediments occur simultaneously) (Figure 1). A diverse range of external forcing, such as tidal pumping, topographically-induced flow, wave pumping and velocity shear, can contribute to the transport of particles across the SWI (Santos et al., 2012; Huettel et al., 2014; Voermans et al., 2018). The chemical concentration gradient is the driving force for molecular diffusion across the SWI (Liu et al., 2014). Small-scale transport processes at the SWI include molecular diffusion, dispersion and turbulent mixing (Santos et al., 2012).At present, six methods are used to study mass exchange flux at the SWI: the concentration diffusion model of interstitial water, the static culture of the original column, the flow culture of the primary column, in-situ box observation, the mass conservation model and flume experiments. Among these, the static culture of the primary column, the low culture of the original column and the in-situ box observation mostly account for the consistency in experimental conditions and field environmental conditions. However, in-situ box observation cannot be used to conduct long-term scaled experiments, and this is primarily because oxidation-reduction reactions will significantly change in closed-loop environments. The mass conservation model is mostly used in closed-looped water or simulated water environments, which belongs to the “black box” model. This method is difficult to employ because it needs to accurately measure changes in nutrient inputs and outputs in a specific volume of water. Therefore, the in-situ box observation and the mass conservation model method are seldom used in the calculation of substance flux. Below is a brief introduction of the four commonly used methods2.4.1 Static culture of the original columnThe intact core sediment incubation device developed by the China Institute of Water Resources and Hydropower Research (IWHR) can be used to collect original column samples. For each sample, sediment samples at a 0-20 cm depth and overlying water at a 0-20 cm depth are collected. Following sample collection, the original column sample collected in the field should be transported back to the laboratory, while the overlying sample of the column sample is carefully removed through siphoning. After algae and suspended material are removed by filtration, the water sample is slowly added to the original sediment column along the pipe wall. For each sample, a syringe is used to extract 50 mL of the water sample at a depth of 20 cm and to measure the elemental concentration in the water sample. After measuring the water sample, the same amount of the overlying filtered water sample is slowly added to the water sample along the tube wall (Wang et al., 2018). The release rate is calculated as follows:

where F is the release flux (mg/m²/d); V is the volume of overlying water in the column (L); c_n, c₀ and c_j-1 are the concentrations of a particle (mg/L) at the nth, 0th and (j-1)th sampling point; c_a is the concentration of a substance in the amended (added) water sample (mg/L); V_j-1 is the volume of the (j-1)th sample (L); S is the contact area of water-sediment in the column (m²) (Figure 3); t is the release time (d). A positive release flux indicates that material within the sediment is released upwards, and a negative value indicates that material within the overlying water is adsorbed by the sediment (Cowan and Boynton, 1996). The biggest disadvantage of the static culture of the original column is that the release or adsorption of material in overlying water cannot be constantly maintained. In addition, this method is also constrained by the sidewall effect of the sampling tube (Xu et al., 2009).

2.4.2 Concentration diffusion model of interstitial waterThe concentration diffusion model of interstitial water is the most commonly used method for calculating material flux at the SWI. If the diffusion of interstitial water obeys the first order reaction of kinetics, the change in its content with depth should obey the law of exponential distribution. Many relevant methods exist, which have developed alongside passive sampling and high resolution sampling technology, to determine the concentration of stratified substances in sediment and overlying water (Liu et al., 2014). However, the peeper device is more accurate and precise. This device is made of organic glass and comprises of a series of chambers (Figure 4). Each chamber side is covered by a dialysis tubing membrane that is biologically inert and has a pore size of 0.45 mm. Ultrapure water is placed into these chambers prior to analysis. Equilibrium of soluble ions and molecules on either side of the dialysis tubing membrane is reached after diffusion through the membrane. The volume of each chamber and the space between chambers can be changed according to experimental requirements. Peeper devices can obtain the complete interface profile between overlying water and pore water. However, the disadvantage of the peeper device is that samplers must be set in sediments for a certain period of time to reach equilibrium before pore water can be sampled (Yang et al., 2016).

The substance concentration of the disturbed layer (3-4 cm) and the overlying water (3 cm) is fitted exponentially at a corresponding depth, and then Fick's first law is used to calculate molecular diffusion flux:

where F is the sediment-water interface diffusion flux (mg/m²/d); $\emptyset $ is the porosity (%) of sediment; D is the effective molecular diffusion coefficient (m²/s); $\frac{{\partial C}}{{\partial Z}}$ is the substance concentration gradient between interstitial water and overlying water (mg/L/cm) (Li and Gregory, 1974; House and Denison, 2002; Zhang et al., 2006; Niu et al., 2013; Yang et al., 2016; Zhao et al., 2017; Wang et al., 2018; Zhao et al., 2018). The concentration diffusion model of interstitial water requires a highly precise interstitial water concentration gradient and is greatly influenced by biological disturbances and sediment stability (Xu et al., 2009).

2.4.3 Flow culture of the primary columnFor this method, a piston with a rubber ring should be inserted at the upper end of the organic glass tube to form an air-tight environment at the upper end of the sediment column. The overlying water is then fully mixed by circulating inlet and outlet water together. The interface diffusion flux formula is as follows (Rydin, 2000; Xu et al., 2009):

where F_n is the diffusion flux of the nth sampling (mg/m²/d); c_n and c₀ are the concentrations of a particle (mg/L) at the nth and 0th sampling point; V is the peristaltic pump flow rate (mL/min); S is the SWI boundary in the column (m²); 60 and 24 are the time conversion factors (Figure 5). Through constant water flow, the flow culture of the primary column method can maintain nutrient content, dissolved oxygen content, redox condition, acidity, basicity and the temperature of sediments and overlying water at a constant level throughout the experimental stage, which safeguards environmental conditions of the culture system when sampling.

2.4.4 Flume experimentThe wave flume test analyzes variation characteristics in substance concentrations under hydrodynamic disturbances, thus directly reflecting the influence of disturbances on sediment- water interface diffusion flux. Tap water is used with sediment to balance overlying water in the flume over a period of several days. One end of the tank is equipped with a push- plate wave maker, generating waves at a height of 3-22 cm and for a period of 0.8-1.5 s (Figure 6). Prior to the experiment, sediment must be disturbed sufficiently and then allowed to settle naturally. When the bulk density of the sediment in the 2 cm surface layer reaches 1.32 g/cm³, the wave maker should be turned on (Zhu et al., 2005). Compared to the static culture test, the flume test uses a larger test system, which can better resolve the sidewall effect, although using undisturbed sediment as a test material is difficult (Xu et al., 2009; Gou et al., 2020).

In recent years, advancements in monitoring methods have improved the accuracy of monitored substances. Moreover, the construction of field stations has greatly increased both the time density and spatial density of ecological monitoring. Rapid improvements in computer performance has also increased simulation accuracy of large-scale distributed hydrological basin models. This coincides with the increasing attention that researchers are paying to ecological environments, the increased demand for flux data as well as the strict accuracy and precision requirements of data. These factors make the selection of the calculation method progressively more important. However, the monitoring and calculation of mass flux has not been standardized or streamlined, and conventional and new calculation methods are being used simultaneously.Therefore, to achieve a balance between efficiency and accuracy and to reach compromises between the standardized methods and the application of new technology, before we estimate the flux, we should fully consider the experimental conditions, the use of the estimated results and other practical situations. The relationship between these different methods must be clarified while also taking into account the advantages and disadvantages of different estimation methods (including accuracy, precision, workload, cost and other factors) (Table 2) in choosing the most suitable method for research.