The Amur River
The Amur (or Heilong Jiang in Chinese) is the world’s 10th largest river with a total drainage basin of ca. 1.89 million km2 and a total length of ca. 4440 km. The vast majority of the basin is located in Russia (53%) and China (45%). Mongolia hosts the remaining 2% of the basin area (Fig. 1). Over a total distance of ca. 2500 km, the Amur River forms the border between China and Russia, and on this entire stretch of the river, no in-situ discharge observations are available. The Amur River ultimately drains into the Tatar Strait between the Sea of Okhotsk and the Sea of Japan. The Amur River is a global biodiversity hotspot hosting endemic fish species and large migratory fish populations as well as huge wetland systems (Egidarev et al., 2016; Simonov et al., 2019). While floodplains on the Chinese side of the river have been severely affected by river regulation (Jia et al., 2020), wetlands in the Russian portions of the basin remain largely intact.
Because it is located in a latitude range from 41 to 56 degrees north, the basin is dominated by cold continental weather with dominant snowfall in winter. Large portions of the river are ice-covered during the winter months. Ice cover monitoring using satellite imagery and satellite altimetry datasets (Zakharova et al., 2021) confirms that the river is frozen from end-November to end-April. The Amur River has several large tributaries (Fig. 1), the largest of which is the Songhua River, joining the Amur from the right-hand side near the town of Tongjiang in China. There are 19 large dams in the Amur River Basin (Simonov et al., 2019). Seven reservoirs have a storage capacity larger than 1 km3, of which two are located in Russia and five in China (Fig. 1). Flooding is common in the Amur basin and seasonal and inter-annual variations of river flow can be related to large-scale atmospheric patterns (Tachibana et al., 2008). The most recent disastrous flood occurred in 2013 (Danilov-Danilyan et al., 2014) and another large flood occurred in 2019. The evolution of flood risk in a changing climate is of concern (Nohara et al., 2006; Yu et al., 2013).
Rainfall-Runoff Modelling
In order to estimate spatio-termporally distributed runoff forcings for the hydraulic model of the Amur River, we set up and calibrate a basin-scale rainfall-runoff model, because available in-situ discharge records are sparse and unevenly distributed. Kalugin and Motovilov, 2018 report the only basin-scale rainfall-runoff modelling effort for the Amur in the open literature. We used the NAM rainfall runoff model (Nielson and Hansen, 1973), which is integrated into DHI’s Mike Hydro River package, for rainfall-runoff simulation. The Amur River basin was divided into 43 individual subcatchments (Fig. 1), using the MERIT hydro DEM (Yamazaki et al., 2019) and the hydrographic DEM processing software TauDEM (Tesfa et al., 2011). The rainfall-runoff model does not include the areas contributing to Lake Hulun in Mongolia, which is essentially endorheic and only occasionally overflows into the Argun River.
As precipitation forcing for the NAM rainfall-runoff model, we used NASA’s Global Precipitation Measurement (GPM) Integrated Multi-satellitE Retrievals for GPM (IMERG) product, more specifically the final precipitation L3 half hourly 0.1 degree x 0.1 degree product, version 06 (Huffman et al., 2019), aggregated to daily values. IMERG precipitation was evaluated against a few available in-situ precipitation stations using straightforward grid-to-point comparison. Resulting double-mass plots are shown in Fig. 2. Although the double-mass plots indicate inconsistencies and shifting biases between the stations and the IMERG product, we concluded from this analysis and from the published literature (e.g. Jiang and Bauer-Gottwein, 2019) that the IMERG precipitation product is suitable for large-scale hydrologic modelling in this region. Gridded land surface (2m) temperature estimates were obtained from ERA5-Land hourly data, provided through the Copernicus Climate Data Store (Muñoz Sabater, 2019). Hourly temperature data were aggregated to daily maximum, minimum and average temperatures. Daily temperature statistics were used to estimate reference ET using the approach by Hargreaves and Samani, 1985. In the NAM model, daily average temperature further controls snow accumulation and snow melt via a threshold temperature for snow fall and a temperature index parameterization of snowmelt (Hock, 2003).
For the 10 in-situ river discharge stations reported in Table 1, the NAM model was automatically calibrated assuming uniform parameters across the entire subcatchment corresponding to the station. In total, 9 NAM parameters (Umax, Lmax, CQOF, CKIF, CK1,2, TOF, TIF, TG and CKBF, please refer to Nielson and Hansen, 1973 and Madsen, 2000 for a description of NAM parameters) were automatically adjusted between reasonable a-priori bounds to minimize overall root mean square error between simulated and observed runoff and overall water balance error, using a global search algorithm as described in Madsen, 2000. Performance was benchmarked against the mean of all observations using the Nash-Sutcliffe Efficiency (NSE). NSE produces optimistic skill scores for seasonal rivers and we therefore also report a skill score in which runoff climatology (i.e. the average of all historical observations for a given day of the year) was used as the benchmark. The climatology index (CI) is calculated as (see also Bennett et al., 2013)
Table 1
Rainfall-runoff model calibration and validation results
Station
|
Catchment ID
|
IMERG Runoff coefficient
|
Calibration Period
|
Validation Period
|
RMSE (m3/s) Calibration
|
RMSE (m3/s) Validation
|
WBE (%) Cal.
|
WBE (%) Val.
|
NSE Cal.
|
NSE Val.
|
Calibration Climatology index
|
Validation Climatology index
|
Novomikhai-lovka
|
23
|
0.39
|
2008–2014
|
2015–2018
|
37.30
|
66.80
|
0.07
|
-3.70
|
0.57
|
0.49
|
0.17
|
0.42
|
Tynda
|
28
|
0.49
|
2008–2014
|
2015–2018
|
54.30
|
45.90
|
-1.00
|
2.71
|
0.54
|
0.55
|
0.30
|
0.34
|
Zvenievoy
|
22
|
0.45
|
2008–2014
|
2015–2018
|
159.60
|
183.10
|
-0.11
|
-3.70
|
0.65
|
0.68
|
0.31
|
0.49
|
Khor
|
37
|
0.71
|
2008–2014
|
2015–2018
|
341.80
|
351.30
|
22.30
|
20.70
|
0.54
|
0.48
|
-0.05
|
-0.07
|
Gouda
|
24
|
0.63
|
2008–2014
|
2015–2018
|
350.60
|
295.80
|
3.10
|
1.20
|
0.72
|
0.75
|
0.22
|
0.32
|
Ust-Niman
|
26
|
0.45
|
2008, 2010, 2011
|
2012, 2013
|
353.60
|
497.30
|
-2.20
|
17.50
|
0.42
|
0.42
|
-0.15
|
0.33
|
Birobidjan
|
25
|
0.63
|
2008–2014
|
2015–2018
|
98.10
|
97.50
|
11.20
|
7.70
|
0.60
|
0.45
|
0.21
|
0.06
|
Ust-Ulma
|
27
|
0.49
|
2008–2014
|
2015–2018
|
776.60
|
771.10
|
2.43
|
-7.00
|
0.52
|
0.48
|
0.01
|
0.17
|
Krasnoya-rovo
|
15
|
0.26
|
2008–2014
|
2015–2018
|
104.60
|
101.80
|
27.60
|
3.30
|
0.51
|
0.76
|
0.08
|
0.69
|
Dalai
|
33 + 34
|
0.12
|
2016–2018
|
2019
|
183.60
|
427.70
|
0.29
|
3.63
|
0.80
|
0.85
|
0.62
|
0.83
|
Mean
|
|
|
|
|
|
|
|
|
0.59
|
0.59
|
0.17
|
0.36
|
where RMSENAM is the root mean squared error between the observations and the NAM simulation and RMSEClim is the root mean squared error between the observations and the runoff climatology.
Transfer of NAM parameters to ungauged subcatchments was based on catchment similarity, using an approach described in Kittel et al., 2020. We used average rainfall, average temperature and average terrain slope as the attributes defining similarity. Ungauged catchments inherited parameters from the gauged catchment that was closest in terms of total normalized distance between the attributes of the two catchments. The standard deviations of the attributes across all 43 subcatchments were used to normalize the distances. Parameter transfer relationships between catchments are illustrated in Fig. 3.
Processing of ICESat-2 land elevation datasets
ICESat-2 is a spaceborne green lidar mission (532nm), mapping the Earth’s surface at unprecedented spatial resolution of approx. 70 cm along track since 2018 (Markus et al., 2017). ICESat-2 is configured with 3 beam pairs that allow for across-track slope determination (90m between pair members and 3.3km between pairs). Each beam pair includes a strong beam (right with respect to orbit direction) and a weak beam (left with respect to orbit direction) with a power ratio of 4:1. ICESat-2 is on a 91-day repeat orbit, for this reason the ground sampling pattern is very dense, while the temporal resolution is low. We used two different ICESat-2 data products: ATL08, which is a low-resolution terrain elevation product (Neuenschwander and Pitts, 2019) and ATL03, which is the native-resolution geolocated photon product (Neumann et al., 2019). We used version 5 of both products and accessed the data through the online portal of the US National Snow and Ice Data Center (https://nsidc.org/data/atl08/versions/5).
The Amur riverbed geometry was extracted from ICESat-2 elevation transects across the Amur River. ATL08 ground tracks were manually inspected to find crossings, spaced approximately every 10–20 river-km, with sufficient data density over the area of interest, which were directed approximately perpendicularly to the river centerline. The ICESat-2 laser is sensitive to cloud cover and mist, which strongly reduces the number of crossings that can be used to extract cross-section geometry. Moreover, we only used crossings during the dry season (i.e. November to March), when water levels in the river are low and a large portion of the cross section geometry is exposed and observable by ICESat-2. The selected ATL08 transects were used to pre-filter the corresponding ATL03 transects and ATL03 data points with elevations outside ± 10m of the interpolated ATL08 elevation were rejected. Remaining ATL03 data was smoothed with a Savitzky-Golay filter (Savitzky and Golay, 1964) using a 3rd degree polynomial to fit the ATL03 points and a variable window length to achieve appropriate smoothing of the ATL03 point cloud. Variable window length was necessary, because cross sections had different absolute length, ranging from hundreds of meters in the upstream portions of the river to tens of kilometres in the large floodplains, and because ATL03 point density varied greatly with atmospheric conditions. All ATL03 points falling more than 1m from the filtered line were removed. The surface elevation geometry was created with a smoothing spline function from the remaining ATL03 points. The degree of smoothing was controlled manually for each cross section to achieve an appropriate representation of the elevation profile. Using the spline interpolation, cross section geometry was resampled to 5m spatial resolution. Open water surfaces were identified in the cross section as entirely flat and smooth surfaces. In the ATL03 datasets for the Amur cross sections used here, we have been unable to detect useable returns from the submerged riverbed. In some sections, we see scattered photons returned from below the water surface, which may be reflected from the riverbed, but the signal-to-noise ratio is too low to enable robust retrieval of submerged riverbed geometry. For this reason, submerged riverbed elevation was extrapolated using the power channel model (Leopold and Maddock Jr., 1953). The power channel model can be written as
$$A=\alpha \bullet {d}^{\beta }$$
$$w=\frac{\partial A}{\partial d}=\alpha \bullet \beta \bullet {d}^{\beta -1}$$
Equation 2
where A is the flow cross sectional area, w is the flow surface width, and d is the flow depth. The parameters α and β are empirical fitting parameters. We assumed a uniform value of the shape parameter β (= 0.2). Depths were estimated as 0.7 times the bankfull depths reported in Andreadis et al., 2013 and available online at http://gaia.geosci.unc.edu/rivers/. In the downstream reaches of the Amur River, which are affected by backwater from the ocean, depths for ICESat-2 cross section acquisition dates were assumed to be equal to the bankfull depths reported in Andreadis et al., 2013. Parameters α were subsequently determined for each cross section from the assumed depth and β and the observed flow width from ICESat-2. Once the parameters of the power channel model were determined, we estimated the submerged riverbed elevation at 5m along-track spacing from the power channel model, and prepared the final cross section for input into the hydraulic model. This included tagging each cross section with the corresponding river chainage, the angle of intersection with the river centerline and sorting the elevations in the direction from left bank to right bank. For selected cross sections, a priori estimates of depth were subsequently manually updated to match simulated spatio-temporally distributed water surface elevations to observed WSE data from in-situ stations and satellite radar altimetry, see section on model validation below.
Hydraulic modelling
Hydraulics in the main branches of the Amur and Songhua rivers (thick blue lines in Fig. 1) were simulated using the fully dynamic version of the 1-dimensional De Saint-Venant equations. Tributary flow (i.e. thin blue reaches in Fig. 1) was simulated using Muskingum routing (Chow, 1988), assuming a kinematic wave speed of 100 km per day and X = 0.25. The 7 major reservoirs in the basin (Fig. 1) were implemented as storage nodes in the Muskingum routing scheme. Evaporation from the open reservoir water surface was neglected and, in the absence of information on reservoir operation, the regulated outflow was determined using a standard operation policy (SOP, Maass et al., 1962) with the target release equal to the annual average runoff and the flood control volume equal to the reservoir storage capacity. This very approximate representation of reservoir regulation will cause significant errors in simulated river flow locally, but, because only a small fraction of total runoff is regulated, the impact on simulated flows in the main Amur and Songhua River reaches is expected to be moderate.
A numerical hydrodynamic model for the main Amur and Songhua River reaches was implemented in the Mike Hydro River software (Havnø et al., 1995), which uses a 6-point finite difference scheme on a staggered grid to solve the coupled continuity and momentum equations (Abbott and Ionescu, 1967). We used a maximum grid spacing of 5 km and a fixed time step of 5 minutes in the numerical solution. The hydrodynamic model was forced with boundary runoff from the rainfall-runoff model and the tributary routing reaches. At the ocean boundary, a constant water level at 0 mamsl was assumed. This will introduce significant model errors locally, because the coastal water level in the Tatar Strait, into which the Amur River flows, is subject to significant tidal variations. However, boundary errors only affect simulated water levels a few tens of km upstream of the boundary. Cross section geometry was imported from the ICESat-2 processing results described in the previous section. We parameterized the friction between the flow and the river bed using Manning’s equation, which expresses the friction slope as dependent on the roughness parameter (Manning’s n), the cross section geometry, and the water level (Chow, 1988). We assumed a global uniform value of Manning’s n equal to 0.033 s/m1/3 during the unfrozen period, except for the most downstream 200-km section of the Amur River, where n = 0.013 s/m1/3 was assumed. Moreover, Manning’s n during the frozen period of the river (end November to end April) was assumed to be 3 times as high as during the unfrozen period and a transition period of 15 days was assumed between frozen and unfrozen states, over which Manning’s n was assumed to vary linearly in time. The factor of 3 between the Manning numbers for frozen and unfrozen states was derived from the inspection of in-situ rating curves prepared by the Russian Hydrometeorological Service, who use different rating curves in the frozen and unfrozen periods (Kouraev et al., 2004).
Validation with satellite-derived and in-situ water surface elevation datasets
In order to validate the hydraulic model and demonstrate its value for water level prediction, simulated water levels were compared with in-situ station datasets from 12 stations (7 in Russia and 5 in China, see Fig. 4) and dozens of satellite altimetry virtual stations (VS, Fig. 4). We included all satellite altimetry time series available on the Hydroweb database (Crétaux et al., 2011, https://hydroweb.theia-land.fr/), for the simulated domain, i.e. 116 virtual station time series in total. For all in-situ and virtual stations, water level time series were extracted from the hydrodynamic model results and were directly compared with the observations. Because we referenced the ICESat-2 cross sections to the EGM2008 geoid model (Pavlis et al., 2012), simulated water surface elevation was also referenced to EGM2008, as were all the satellite altimetry observations at VS. The vertical reference of the in-situ stations was unknown and we therefore expect time-constant bias between the model and the in-situ observations.