Flow dynamics in large tidal delta of the Northern Dvina River: 2D simulation

A tide-dominated river delta with anastomosed distributary network and a vast floodplain area is believed to be one of the most complicated hydrological objects from the point of view of flow dynamics. The high temporal and spatial variability of its characteristics, unsteady flow and reversing currents, secondary turbulent structures and irregular mean flow velocity distribution make tidal delta modelling a challenging task (Mikhailov 1997; Hoitink & Jay 2016). Nevertheless, the fact that many large cities and ports are situated in estuaries and delta areas demands an urgent supply of current data and forecasts of water regime parameters. This is of great importance for navigation, as well as for protection from flooding resulting from storm surges and snow-melt floods.

In major river deltas direct field observations still remain an inadequate research tool despite the fact that new methods, such as satellite positioning, acoustic profiling and automated water level logging, are being actively devised nowadays. These methods are costly and labour-intensive. They require work safety under adverse environmental conditions, which is of great interest to hydrology. However, the data provided by field measurements and observations offer an excellent background for development, calibration and validation of numerical models, which in turn allows simulation of any hydrological situation including the most dangerous and exotic ones.

In the last decade numerical models reproducing hydrodynamics, sediment transport and water mixing were developed for a number of large tidal estuaries: the Delaware estuary on the Atlantic coast of the United States (Aristizábal & Chant 2013), the Knik Arm and the Cook Inlet in south-central Alaska (McAlpin et al. 2013), the Ems estuary between North Germany and the Netherlands (Pein et al. 2014). In these cases, the river runoff was less important for the results of modelling.

An interaction between river flow and tides is much more significant for shallow tidal rivers (Buschman et al. 2009; Hoitink & Jay 2016). Buschman et al. (2010) studied how tides change the flow distribution in the idealized delta using the 2D model of channel bifurcation. By means of numerical experiments with length, depth and bed roughness of distributaries the subtidal (tidally averaged) water discharge allocation was investigated. Sassi et al. (2011) provide the test of the idealized scheme of Buschman et al. (2010) under conditions of a real delta. A 2D model of the Mahakam Delta (Indonesia) reveals others factors affecting subtidal discharge distribution, such as channel pattern and junction angles of branches.

Simulation of Hurricane Gustav's consequences for the lower Mississippi River in 2008 was presented by Martyr et al. (2012). This research managed to improve validation results of the 2D model of Bunya et al. (2010) regarding both the river flow variability during the hurricane event, and variations in the bottom stress depending on flow velocity regime.

Pan et al. (2014) developed a coupled 1D-3D model of the Pearl River delta (China) which was calibrated and validated on the basis of an extensive field dataset. However, the validation is provided only for subtidal water levels, velocities and salinity fields, and variations of hydrodynamic parameters during the tidal cycle were not considered. This lack is not present in a 2D model of the Han River large tidal estuary on the Korean peninsula, which was developed by Park & Song (2017) in order to study pollutant concentration fluctuations in the tidal reach. It was found that reversed flows were repeatedly created and destroyed as the tide varied. The change in the flow direction was related to oscillatory variations of the pollution, and the horizontal recirculation zone at the front of the reversed flow became a storage area where the pollutant cloud was trapped.

The hydrodynamic model of the Lower Northern Dvina River, presented in this paper, is the first 2D model of this area. Our current research could be considered not only as one more case study of a real complex tidal delta, but also as an attempt to evaluate the ability of the model to simulate variations of runoff distribution, flow velocities and water level variations during the tidal cycle as well as tidally averaged flow characteristics.

The Northern Dvina is one of the greatest rivers in the north of European Russia, and the oldest Russian seaport of Arkhangelsk is situated at its mouth. The mouth area of the Northern Dvina River includes an extensive delta with anastomosing distributary network and vast lowlands subject to inundation by tidal and surge water level oscillations, as well as to flooding by excessive river runoff (Zotin & Mikhailov 1965). Tidal waves, wind surges and reversed flows progress along the main delta channels toward the delta top and even further upstream. The island territories of the city of Archangelsk suffer from both storm surges and river spring snow-melt flooding. Regular observations of water level by automatic recorders in the delta have been carried out by the Northern Headquarters of Federal Service for Hydrometeorology and Environmental Monitoring since the 1950s (Zotin & Mikhailov 1965). The river discharge measurements in delta distributaries throughout the whole tidal cycle have been organized several times since the 1960s by field teams from the Zubov State Oceanographic Institute and Lomonosov Moscow State University. The first attempt to develop a numerical hydrodynamic model of the Northern Dvina delta was undertaken by scientists of the Russian Arctic and Antarctic Research Institute (Kotrehov & Pavlova 1983), and was a 1D model based on Saint-Venant equations. The model was expected to reveal transformation of tidal and surge waves propagation along the Lower Northern Dvina River. It was shown that the maximum water level during a wind surge is almost the sum of a range of tidal and surge waves. Results of wind surge simulation using that model were rather accurate in terms of water levels; however, the accuracy of flow distribution by delta distributaries was not considered.

The goal of this paper is to demonstrate the advantages of the numerical model based on STREAM_2D software (Belikov et al. 2002, 2015) for flood routing, e.g., to calculate water level elevation and flow parameters, to delineate the flooding zones and the depths of inundation, and to evaluate river runoff and reversed currents for delta channels under the influence of spring and neap tides, storm surges and changeable river runoff. Contemporary data on water levels at state regime gauges and water discharge measurements in delta distributaries during the tidal cycle acoustic Doppler current profiler (ADCP) were used for model calibration and validation. Model sensitivity tests concerning the accuracy of input data and inner parameters were carried out to explain why the simulated results are somewhat inconsistent with the observational information. The scenario of extreme hazardous storm surge flooding that happened in November 2011 on the outskirts of Arkhangelsk was simulated, as well as other typical and extreme hydrological events (Lebedeva et al. 2015).

The results of the model validation and sensitivity tests can be used for improvements in local flood monitoring systems (Alabyan et. al. 2016) and in further model development.

The Northern Dvina River discharges into Dvina Bay on the White Sea and is the second greatest river in the north of the European part of Russia. It drains a catchment of about 357,000 km² and has a mean freshwater input of 3,330 m³/s. The mean annual peak discharge is 21,200 m³/s, with the great majority of runoff occurring during the spring due to snow melting. The mean low-flow river discharge is equal to 1,600 m³/s in summer and 740 m³/s in winter (Lebedeva et al. 2015).

The mouth area of the Northern Dvina River consists of a multi-channel delta which has a shape resembling that of an equilateral triangle with a surface of approximately 900 km² and about 100 km reach upstream affected by tidal variations. The delta is constituted by three main anabranches: Nikol'skiy, Murmanskiy and Korabel'nyy, and two large distributaries: Maymaksa and Kuznechikha. The main waterway is routed via Maymaksa, characterized by a narrow and deep channel (Figure 1).

Sea level oscillations are both tidal and caused by storm surges. The tidal regime of Dvina Bay is of semi-diurnal type with a period of about 12.5 h. The tidal range is 1.25–1.5 m during spring tides decreasing to 0.6–0.8 m for neap tides (Zotin & Mikhailov 1965; Lebedeva et al. 2015). Tidal oscillations are distinguished by slack tides (locally called ‘manikha’): the cessation of water level rise during the flood tide for 2 or 3 hours. The reversing currents during the low-water season propagate upstream over distances of more than 60 km, and tides progress up to 135 km from the marine coast reaching the inflow of the Pinega River, a major right tributary. Throughout the high runoff stage the influence of tides is negligible. A combination of river discharge coming from the catchment and sea level regime at the sea margin is the main control defining the flow dynamics within the mouth area (Zotin & Mikhailov 1965).

The island parts of the city of Archangelsk, located at the top of the delta, and a number of minor settlements suffer from storm surges. The surge floods are most frequent in autumn, when the north and northwest winds blow from Dvina Bay and the central part of the White Sea. The water level rise during the most hazardous surge floods reaches 1.5 m (Lebedeva et al. 2015). The most recent occurred on November 15–16, 2011 and was the highest for the whole period of observations since 1884 at Solombala gauge.

The oldest Russian sea and river port of Arkhangelsk was founded at the top of the Northern Dvina delta in the 16th century. Therefore, the hydrography and hydrological regime of the lower Northern Dvina is the most well-investigated among the estuaries of the Russian Arctic. A special river-mouth hydrometeorological station was constructed in Arkhangelsk in 1955. It controlled the work of more than ten hydrological gauges in the river mouth area. Nowadays, there are seven functioning water level gauges (Figure 2), and six of them are equipped with automatic tide-gauges with measuring frequency less or equal to 15 minutes (gauges of Northern Headquarters of Federal Service for Hydrometeorology and Environmental Monitoring and ‘Rosmorport’ industry). At the hydrological station Ust-Pinega, located just downstream from the junction with the Pinega River, the water level is measured twice a day. The river runoff is measured there several times a year and these measurements are used to update stage-discharge curve. The location for Ust-Pinega station was chosen as the place where tidal effect is usually not significant (Zotin & Mikhailov 1965). However, extreme wind surges can achieve this point (Kotrehov & Pavlova 1983).

Special field studies were carried out in the Northern Dvina delta in the 1960s and 1980s. They included measurement of the water discharge distribution between the main branches and distributaries in different river flow conditions. As a result of data assimilation and analysis the correlation between residual discharges in the main branches and the total river discharge at the top of the mouth area (Ust-Pinega station) was established (Polonskiy & Kuz'mina 1986). Nevertheless, the influence of sea level change and the tidal range impact on the flow distribution amid the branches was not taken into account properly. The importance of these factors becomes evident in the model sensitivity test (see below). During the summer seasons in 2011 and 2013, water discharges were measured at six cross sections in the main delta branches throughout the entire tidal cycle (approximately 12.5 hours) using ADCP (Polonskiy & Mishin 2013; Lebedeva et al. 2015). Measurements were taken in each of the six cross sections (Figure 2) during the full tidal cycle and were not synchronized. Most of these data were used for the model calibration and validation.

STREAM_2D is a software package for two-dimensional numerical simulations. It has been developed over the last three decades by V.V. Belikov and V.V. Kochetkov and was certified by the Russian Federal Service for Intellectual Property in 2014. STREAM_2D is well adapted to resolve problems of flow dynamics in rivers with complicated channel pattern and wide floodplain, even with a lack of field data. It has been tested on a large number of Russian rivers (Belikov et al. 2002, 2015; Alabyan et al. 2014) but has never been used before for tidal estuaries with reversing currents. We chose STREAM_2D in an attempt to assess its suitability to simulate the hydrodynamics of large anastomosed deltas subject to tides.

In contrast to 2D diffusive wave flood models (also called zero-inertia) in which it is possible to simplify SWE by neglecting the inertia terms in the momentum equations (Jahanbazi et al. 2017), the flow pattern in tidal estuaries can be simulated only by the fully dynamic equations because the inertia features are of key importance there. The numerical solution of these non-stationary two-dimensional equations in a STREAM_2D model is carried out using a Godunov-type method based on an exact solution to the Riemann problem for the shallow-water equations (Belikov & Semenov 1997; Aleksyuk & Belikov 2017; Norin et al. 2017). A flexible hybrid triangular–quadrangular computational mesh is applied to cover a modelled area. Quadrangles describe channels and triangles constitute river floodplains and delta lowlands. Such kinds of mesh seem to be the best way to describe topography and bathymetry in a complicated channel network with adjacent submerging terrains.

As the current version of the Northern Dvina mouth model has been developed focusing on the simulations of tidal and synoptic timescale processes, the riverbed deformations and sediment transport by the river and along the seashore, as well as the effects of density variations in the mixing zone, are not taken into account. In fact, the intrusion of salt water into the channel of the Maymaksa distributary and possibly other deeper channels of the delta, occurs regularly at a distance of about 15 km from the marine edge it is possible to form a wedge of sea water with salinity of up to 15‰ in a water layer from the bottom to 5 m above it (Lupachev 1976). In extreme cases, the salt-water intrusion can reach the delta top and disturb the normal operation of Arkhangelsk water-intakes (Magritsky et al. 2017). However, on the contemporary modelling stage the water flows are assumed to be of homogeneous density.

The upper open boundary of the model (‘river’) is located at the Ust-Pinega hydrological station. Tidal oscillations are assumed not to affect the river flow at this distance (Zotin & Mikhailov 1965). The river discharge at Ust-Pinega (and its variation in time for non-stationary cases) was set as the upper boundary condition. The lower boundary (‘sea’) of the model area is established at 2–3 km seaward from the delta marine edge. Time series of water level oscillations at the gauges of Severodvinsk and Mud'yug were used to define the lower boundary condition (Figure 2).

A computational flexible mesh was generated using the mesh-making tools of the STREAM_2D package. The mesh consists of more than 190 thousand cells (Figure 3(a)). The linear dimensions of cells range from 10 m in channels of minor distributaries to 100–200 m over the floodplain and delta coastal lowlands.

The bathymetry of channels was defined by up-to-date nautical charts and hydrographical surveys of waterways and navigational passes. Terrain relief was taken from topographic maps and the mesh was constructed in accordance with the resolution of bathymetry data. The collected relief data were digitized, represented as a field of vector points, and then interpolated to the centres of the mesh cells to construct the digital elevation model (DEM) of the Northern Dvina River mouth area.

The calibration process consisted of searching for the field of roughness coefficients n (x, y), which reproduce a real hydrological situation most accurately. In other words, it is a situation in which modelled water levels at gauges and water discharge distribution between branches have the best correspondence to the measured values.

The model area was divided into ‘roughness zones’. Each zone has a unique roughness coefficient value. The floodplain and the delta surface are united into one zone. The boundaries of this zone are defined based on the vegetation contours visible on Landsat satellite images. Channels of the five main anabranches and distributaries of the delta represent zones of different roughness (Figure 3(b)).

The type of flow dynamic conditions in the Northern Dvina delta significantly varies depending on the river runoff regime. The ‘non-tidal’ period of a year occurs throughout spring snow-melt flood; when river discharge is high the tidal wave does not progress in the upstream direction and the flow does not turn backwards over the entire delta during the whole tidal cycle (hydrology of the Northern Dvina River mouth area 1965). As for the rest of the year, the depth and cross-sectional areas of channels increase during flood tides and the flow direction can be reversed. This type of regime can be defined as ‘tidal’.

A channel reach can be characterized by different roughness under tidal and non-tidal periods. Therefore, the model should be calibrated separately for low-water (tidal) and high-water (non-tidal) conditions. Two sets of roughness coefficients n (x, y) were found as a result of calibration (Table 1).

The model calibration was based on the data of hydrological monitoring and the discharge analysis throughout July 1983 for tidal conditions. The period of peak discharges in April 1983 was considered to calibrate non-tidal conditions. In 1983, the largest number of water level gauges worked simultaneously (five gauges in the delta area and Ust-Pinega station) and discharge measurements were conducted in the delta branches (Polonskiy & Kuz'mina 1986). Roughness coefficient for the river stretch upstream of the delta (from Smolnyy Buyan gauge to Ust-Pinega station) was defined according to the best consistency between the modelled and the observed water levels at Ust-Pinega (Table 1). During the calibration process it was becoming clear that channel roughness mainly affects the water surface slope along the river stretch upstream from the delta and has less effect on determining water level elevation in the delta. Roughness coefficients for the delta channels were selected in order to minimize the error in the calculation of water discharge distribution between branches (Tables 2 and 3). It was not possible to achieve an adequate consistency between the calculated and measured residual distribution of discharges in tidal conditions. The largest error of 17–18% is found in flow distribution between the Korabel'nyy and Nikol'sky branches (Table 2), which could be connected with both the uncertainty in the marine boundary conditions and the insufficiently detailed bottom relief data of these not intensively navigated channels. This challenge will be discussed below in the section devoted to sensitivity tests. In non-tidal conditions, flow distribution between all delta branches is simulated with errors of less than 2% of the total discharge (Table 3), which looks appropriate.

As a result of calibration, the root-mean-square errors (RMSE) of calculated water levels at gauges are concentrated mainly between 5 and 15 cm for both tidal (low-water) and non-tidal (flood) conditions.

Validation of the model for tidal and non-tidal conditions was carried out using some new data. Model validation for tidal conditions was undertaken for the low water period of July 2013 (Figure 4(a)). The modelled fields of water level were compared to the monitoring data from the following gauges: Smolnyi Buyan, Solombala and Ekonomiya. The model adequately reproduced the tidal range as well as the time of high water and low water. RMSE of the calculated water level was equal to 0.05 m for Smolnyi Buyan, 0.05 m for Solombala (Figure 5) and 0.18 m for Ekonomiya.

Within the period from 26 to 31 July 2013 water discharge time series were measured with ADCP during the tidal cycle in five main delta branches and distributaries (Figure 6). Simulation of water discharge change during the tidal cycle was discovered to be most representative for cross sections at the top of the delta (s1, Smolnyi Buyan), as well as for the Korabel'nyy branch (s4) and the upper reach of the Maymaksa distributary (s5). The range of water discharge magnitudes is significantly underestimated by the model for all cross sections except Nikol'skiy branch (s2), where the simulated reverse water flow was higher than the measured one. Therefore, the simulated residual discharge for Nikol'skiy branch is of an even smaller negative value (Table 4). Nevertheless, the main features of the discharge time series, such as a change of flow direction, slack tide, the time moments of peak flood and ebb discharges, were reproduced adequately. Considering calibration results (Table 2) and the results of non-tidal condition modelling (Tables 3 and 5), it seems that the lower the total river discharge, the more errors there are in simulated residual discharges of the delta branches. Further studies of this challenge are provided in the section ‘Sensitivity tests’ below.

Model validation for non-tidal conditions was carried out based on the data of the snow-melt flood period from 1–19 May 2012 (Figure 4(b)). During the first 4 days of this period, ice phenomena were still observed, so simulated water levels were compared to the observed data for the period from 5–19 May, when the river was completely free of ice (Figure 7). RMSE for calculated mean daily water levels at the gauge Ust-Pinega was 0.15 m; for calculated levels of high and low water at Smolnyi Buyan and Solombala it was 0.07 m and 0.12 m, respectively. These values were slightly exceeded and reached 0.3 m and 0.2 m, respectively, for the period of wind surge of 12–13 May when the northwest wind speed was 8–10 m/s with gusts up to 16–20 m/s.

The comparison of the modelled runoff distribution amidst the delta channels for non-tidal conditions with the data of Polonskiy & Kuz'mina (1986) demonstrates appropriate agreement with accuracy within 2% of the total flow (Table 5).

Both inner model parameters and input data contain uncertainty and measurement errors. Special numerical experiments were carried out to evaluate factors impacting model results. The sensitivity stimulations tested the impact of change in bathymetry, marine boundary water levels, roughness coefficients and mesh refinement on modelled water levels and especially on flow distribution between branches. Specification of these factors can improve model performance in future.

Increasing roughness coefficient in all delta channels by 0.002 leads to an increase in the simulated water levels by 0.1–0.15 m and a decrease in the range of tidal fluctuations by 0.05–0.1 m at Solombala and Smolnyi Buyan gauges. The closer to the model marine boundary, the smaller the impact of roughness on the water level elevation.

The change in roughness coefficient for the floodplain from 0.06 to 0.08 (roughness zone No. 9, Table 1) hardly affects model results in flood conditions. Only upstream from the delta along the reach from Smolnyy Buyan to Ust-Pinega does the water level rise to 0.1 m at Ust-Pinega gauge.

The varying roughness of delta branches (Table 1) led to a flow redistribution between them of no more than 8% for low water stage and 3% for flood conditions. Thus, the most important control of both the discharge distribution in the delta and water levels is not the roughness coefficient.

There is high uncertainty in setting boundary conditions along the marine boundary. Mean water level at Mud'yug gauge is 0.2 m lower than in Severodvinsk. In addition, the high and low water at Mud'yug occur about 40 minutes later than in Severodvinsk. The distance between these gauges is approximately 50 km.

Two methods of setting the boundary condition were tested:

• Option 1: water levels at Severodvinsk and Mud'yug gauges were taken as endpoints of the marine boundary line, and the level along it was calculated by linear interpolation between them.

• Option 2: a single value of water level was assumed for the whole lower boundary and was averaged by the synchronously measured values at Mud'yug and Severodvinsk.

As a result of boundary conditions change from option 1 to option 2, the redistribution of residual discharges between the delta branches amounted to almost 35% of the total river runoff in low water conditions (Table 6). Nikol'skiy and Korabel'nyy branches were mostly affected by the change of boundary condition. The distribution of residual discharge between them is much closer to the observed data under option 2 (Table 6). However, a comparison of simulated time series of water discharges in delta branches with the observed data reveals that a single water level at the marine boundary (option 2) results in significant shifts in water discharge oscillations in time. As a result, the time moments of peak discharges and changes in the flow direction are less consistent with ADCP measurements than under option 1. Moreover, setting up a single water level along the entire marine boundary (option 2) intensifies the slack tide and significantly increases the calculated water levels in the eastern part of the delta: by 20−40 cm for flood and by 5−10 cm for low water stage. Thus, the water level at the marine boundary plays a key role in the distribution of water discharges between the delta branches, especially for low water conditions.

As a result, it was concluded that the parameters of tidal oscillations (A and H₀) have a great impact on the distribution of water flow between the delta branches when Q is less than 15,000 m³/s, A and H₀ being of key importance when Q is less than 5,000 m³/s (Figure 8).

Two scenarios were simulated in order to analyse the sensitivity of the model to errors in bathymetry. Bottom elevations in Nikol'skiy branch were decreased by 1 m (option 1) and increased by 1 m (option 2) as compared to the real ones. This aggrading and deepening of the Nikol'skiy branch resulted in an increase and decrease in the magnitude of the tidal fluctuations by 0.1–0.15 m. Moreover, noticeable changes occur in water discharge time series at the flow bifurcation into Nikol'skiy and Karabelny anabranches (Figure 9). As a result of the bottom modification, the range of water discharge tidal fluctuations in the Nikol'skiy branch more than doubled (from 2,500 m³/s to 5,700 m³/s). However, no significant changes in the residual discharge distribution between the delta branches were found. Flow redistribution did not exceed 2% of the total river runoff.

For mesh refinement, reducing the cell size in channels by 1.5−2 times did not result in significant changes either in the simulated water levels at key gauges or in the flow distribution in the delta. Thus, the appropriate cell size seems to be optimal for the current stage of modelling.

Concluding this section, it can be said that, in contrast to river sections located beyond tidal impacts causing backwater effect and reverse currents, in the Northern Dvina delta the roughness is a less important parameter for runoff distribution compared with other factors and circumstances. For tidal estuaries maximum attention should be paid to secure the model with detailed data about the water level oscillations along the marine boundary, which is especially important for multichannel deltas with extended sea coast.

An extreme wind surge occurred in the Northern Dvina River mouth on 15–16 November 2011. It was stronger than any surge in the 20th century and the whole period of observations since 1884 at Solombala gauge. The first extreme peak was recorded at 10:00 a.m. on 15 November and the second peak came during the night of 15–16 November. It exceeded the first one by a few centimetres, reaching the 1.97 m mark at Solombala gauge. The storm surge was caused by a powerful Atlantic cyclone that invaded the White Sea and on 14–16 November induced north and northwest winds over Dvina Bay with speeds from 8 to 13 m/s and gusts up to 23 m/s.

Mean daily river runoff at Ust-Pinega station for this period changes from 2,280 to 2,480 m³/s. Simulated surge water levels were compared to the data from Solombala, Ekonomiya, Smolnyi Buyan gauges and the Ust-Pinega station. The difference in values of maximum surge water levels obtained through the model and observations are no higher than 0.12 m for all gauges (Figure 10). Time moments of maximum water level at Ekonomiya, Solombala and Smolnyy Buyan gauges were reproduced accurately. Unfortunately, the exact time of surge peak at Ust-Pinega was not registered, and only the water level value was fixed and compared to the model results (Figure 10).

The simulated time series of water levels show that the surge must have two peaks of similar height corresponding to high water of the two successive tidal cycles. Despite the similarity of the water level values, the gradients of water surface and flow velocities at these time moments were very different (Figures 11 and 12). During the first two tidal cycles (since 20:00, 14 November) a wind surge was developing: the water was being pumped into the river mouth area. At the moment of high water at the marine boundary (15 November, at 10:00) there was a significant reverse gradient of water surface, providing reversed currents and large negative discharges in the delta branches. The propagation of the second tidal wave was accompanied by much lower current velocities over the whole river mouth area, because the entire storage capacity of the delta channels was filled and water surface was close to horizontal (Figure 11). The negative water discharges in all branches and at the top of the delta (Smolnyy Buyan) had significantly lower magnitude than during the first peak (Figure 12). Afterwards, the wind surge weakened throughout subsequent tidal cycles, reversed currents stopped, and positive (directed to the sea) flow in the delta gradually accelerated (Figure 12). As a result, the entire water mass pumped into the river mouth was discarded into the sea.

During the wind surge the delta terrain was flooded more than usually happens in spring. Reverse gradient of the water surface provided high water levels in the centre of Arkhangelsk and over the island lowlands. The most significant inundation occurred in the vicinity of the port of Ekonomiya and settlements along the Maymaksa and Kuznechikha delta branches (Figure 13). Pontoon ferries were closed. Significant damage to the city of Severodvinsk was caused by flooding of houses, wave impact and strong winds.

The flood area contours the maximum inundation (Figure 13) and any time of the surge flood can be generated based on the model results. The precision of contours significantly depends on the accuracy of the DEM. Further detailing of the terrain relief can be the next step in obtaining much more accurate modelling contours of the flooding area without changing the modelling procedure itself.

Another way to improve the accuracy of wind surge simulation would be to take into account the speed and direction of wind over the water surface. The current version of the model does not take into account any additional wind effects on the flow. It is assumed that the wind surge wave is being formed completely outside the model boundaries (that is taken into account via low boundary condition), and inside the river mouth area the surge wave transforms due to the interaction with the bottom and coasts. The 1D model of the surge wave propagation in the North Dvina River mouth area (Kotrehov & Pavlova 1983) demonstrated that introducing the additional friction coefficient based on wind parameters could improve the model, but not in every case, because there is not enough adequate data for calibration of the wind friction coefficients.

As observed in the previous sections, the presented model of the Northern Dvina River mouth area simulates water level fluctuations quite accurately. RMS error is less than 0.18 m, which is comparable with the RMS errors of subtidal water level calculations by similar models of the Lower Mississippi (Martyr et al. 2012), the Pearl River delta (Pan et al. 2014) and the Mahakam delta (Sassi et al. 2011).

The developed model precisely calculates the flow distribution between the delta branches in the case of high river runoff and non-tidal conditions. For the validation scenarios, the difference between calculated and observed discharge distribution is less than 2% of the total. However, the errors increase significantly for the periods of low flow in tidal conditions. A similar effect was analysed by Buschman et al. (2010) and Sassi et al. (2011). A number of scenarios and numerical experiments were provided to confirm that the tide range and tidally averaged water level affect deeply residual discharge distribution. It is supposed that this problem is related to a lack of knowlege about water level elevation and its fluctuations along the lower model boundary. Two methods are proposed for solving this problem for further model development: (1) to monitor water levels along the delta marine edge by additional temporal gauges; or (2) to shift the model marine boundary further into the sea and include the entire Dvina Bay, or part, into the model domain. In the latter case, additional bathymetry data are needed to simulate how the tidal wave would transform. An alternative way to enhance the simulation strategy is to arrange the interaction of the estuary model with some hydrodynamic model for the entire White Sea. In this case, the ‘marine model’ would supply the ‘river model’ with lower boundary conditions.

Other factors affecting the simulation results are the channel bathymetry, the floodplain topography and the roughness. Numerical experiments, consisting of varying the bottom elevation in the Nikol'skiy channel, demonstrated the change of the tidal discharge range, but not the residual discharge distribution. Sensitivity tests of the roughness coefficient showed that the less the river runoff and the closer the delta marine edge, the lower the results' sensitivity to Manning's value.

Moreover, the Manning coefficients for the delta branches under the tidal conditions seem to be extremely low. This phenomenon could either be due to numerical scheme specifics or have the physically based nature connected with features of the turbulent viscosity in reversing currents (Alabyan & Panchenko 2017), which could be the subject of special research. In any case, significance of the roughness and the bathymetry of delta arms appeared to be much less important than the water level regime along the marine boundary.

Previously, the residual discharge distribution in the Northern Dvina delta was considered to be a function of the total river runoff only (Zotin & Mikhailov 1965; Polonskiy & Kuz'mina 1986). However, our research shows that it is necessary also to take into account the tide parameters at the delta marine edge, as well as wind surges.

The presented model allows users to simulate floods caused by both wind surges and by spring snow-melting for the total area of the vast delta and floodplain of the Northern Dvina River. More detailed land surface elevation data could enhance the model in terms of more accurately reproducing the inundated zone boundary and the water depths over flooded terrain. However, it should not significantly affect the precision of water level simulation, because the main flow concentrates in channels.

As a result, the STREAM_2D modelling software appeared to be a suitable instrument to simulate the hydrodynamic regime of large anastomosed tidal deltas. Integration of such computer models with field observations can represent a very productive approach for scientific research into flow dynamics in tidal estuaries and also for the practical issues, such as monitoring, forecasting and estimation of hazards, providing hydrological information for navigation, industrial and urban water management, etc.

We are going to incorporate the recently developed model of the Northern Dvina River mouth with intelligent information systems for operational river flood forecast (Alabyan et al. 2016). Therefore, the model can be used to solve water management problems at the Northern Dvina delta.

The research was supported by the Russian Foundation for Basic Research (Project No. 16-05-01018 ‘Research of dynamics of tidal waves and surges in river mouths of the White Sea basin’) in the field data collection for calibration and validation of the model, and the Russian Science Foundation (Project No. 14-17-00155 ‘River runoff parametrization for identification of hydrological hazards and their environmental consequences’) in the model sensitivity tests and wind surge simulation.