Integrated Hydrologic-Hydrodynamic Inundation Modeling in a Groundwater Dependent Tropical Floodplain

The rapid development of free and open-access hydrological models and coupling framework tools continues to present more opportunities for coupled model development for improved assessment of floodplain hydrology. In this study, we set up an Upper Zambezi hydrological model and a fully spatially hydrological-hydrodynamic coupled model for the Barotse Floodplain using GLOFRIM (GLObally applicable computational FRamework for Integrated hydrological– hydrodynamic Modelling). The hydrological and hydrodynamic models used are WFLOW and LISFLOOD-FP, respectively. The simulated flows generated by the wflow model for the upstream gauge stations before the Barotse Floodplain were quite similar and closely matched the observed flow as indicated by the evaluation statistics; Chavuma, nse = 0.738; kge = 0.738; pbias = 2.561 and RSR = 0.511; Watopa, nse = 0.684; kge = 0.816; pbias = 10.577 and RSR = 0.557; and Lukulu, nse = 0.736; kge = 0.795; pbias = 10.437 and RSR = 0.509. However, even though the wflow hydrological model was able to simulate the upstream hydrology very well, the results at the floodplain outlet gauge stations did not quite match the observed monthly flows at Senanga gauge station as indicated by the evaluation statistics: nse = 0.132; kge = 0.509; pbias = 37.740 and RSR = 0.9233. This is mainly because the representation of both floodplain channel hydrodynamics and vertical hydrological processes is necessary to correctly capture floodplain dynamics. Thus, the need for an approach that saves as a basis for developing fully spatially distributed coupled hydrodynamic and hydraulic models’ assessments for groundwater dependent tropical floodplains such as the Barotse floodplain, in closing the gap between hydrology and hydrodynamics in floodplain assessments. A fully coupled model has the potential to be used in implementing adaptive wetland management strategies for water resources allocation, environmental flow (eflows), flood control, land use and climate change impact assessments.


Introduction
Floodplains play an important role in global hydrological and biogeochemical cycles. Many socioeconomic activities depend on water resources in floodplains, which are important components of many large river systems in Africa [1,2]. Modeling floodplain wetland processes is a topical issue in environmental studies. This interest has arguably grown largely due to the need to protect such aquatic ecological systems from the adverse effects of land use and climate change, which are linked to environmentally degenerative anthropogenic activities [3]. In the absence of adequate and continuous historical observations, advances in remote sensing, computing power, and model coupling frameworks have made distributed numerical models an increasingly attractive solution where spatial understanding of floodplain processes is required [4]. For example, using GLOFRIM (GLObally applicable computational FRamework for Integrated hydrological-hydrodynamic Modelling), it is possible to couple models spatially explicit (on a grid-togrid basis) and online (on a timestep basis). This kind of technological advancement enables the user to create a coupled hydrological and hydrodynamic model [5]. For example, by using GLOFRIM, Hoch et al. (2019) [6] coupled the global hydrologic model PCR-GLOBWB (PCRaster Global Water Balance) [7,8] with the hydrodynamic models CaMa-Flood (Catchment-based Macro-Scale Floodplain) [9] and Lisflood-FP [10]. The results of these studies show that replacing the kinematic wave approximation of the hydrologic model with the local inertia equation of CaMa-Flood greatly enhanced the accuracy of peak discharge simulations. Whereas, flood maps obtained with Lisflood-FP improved representation of observed flood, related to downscaled products of PCR-GLOBWB and CaMa-Flood. Results from these studies in the Amazon and Gangas basins confirm that model coupling can indeed be a viable way forward towards more integrated simulations of flooding processes. This is due to the fact that coupling of hydrology and hydrodynamics in inundation models allows for physically more integrated assessments and to compensate for their respective shortcomings [5].
In recent years, the demand for an understanding of the hydrological and hydrodynamic processes for the Barotse floodplains is ever increasing, especially with the advent of climate change and variability and its potential impact on hydrological variables [11][12][13][14]. Spatio-temporal discharge variability and basin characteristics play a significant role in the inundation extent of the Barotse Floodplain [15] and, thus, they control habitat conditions of the floodplain river channels and the linked floodplain, which impacts ecosystem services. Additionally, groundwater in the river channelfloodplain exchange may affect the volume and inundation dynamics in the Barotse Floodplain [16]. However, these interactions are poorly understood despite the fact that they are crucial to understanding the impacts of future changes in the upstream basin and the wetland itself on water availability, predictability of low flows and floods, and possible threats to the ecosystem and the unique socio-economical system supported by it. Currently, for the Barotse floodplain, modelling approaches to try and understand the floodplain dynamic processes are based on one-way coupling [2,16,17]. However, as Hoch et al. (2017b) [18] reason, the resulting lack of interaction between hydrology and hydrodynamics, for instance, by employing groundwater re-infiltration on inundated floodplains, can hamper modelled inundation and discharge results where such interactions may potentially be important. In this study, we set up a hydrological model for the upper Zambezi Basin and a framework for the development of a fully spatially coupled hydrological-hydrodynamic model for the Barotse Floodplain. This is significant in supporting the development of coupled models in order to understand the multi-way interactions between river flows, wetland inundation, and groundwater as these processes still remain poorly understood in the Barotse Floodplain.

Study Area
The Barotse Floodplain is found within the Upper Zambezi Basin (UZB). The Upper Zambezi is the broad extent of the Zambezi River from the source 25 km southeast of Kalene Hill in Mwinilunga District, North-Western Zambia, through Angola and Barotse Floodplain to the Victoria Falls. The basin lies between latitudes 11°S and 19°S, and longitudes 18°E and 27°E, which covers part of western Zambia. While the Barotse Floodplain lies between 14°S and 17°S, and longitude 22°E and 24°E ( Figure 1). The flooding in the Barotse floodplain is a consequence of the breaking of the banks of the Zambezi River, reaching the peak in April. The floodplain measures approximately 240 km long, stretching from Lukulu to Senanga and 34 km wide. The total floodplain area is estimated at 7,700 km 2 [19]. The floodplain is a sanctuary to variety of biodiversity both flora and fauna [15,20,21].

Methodology
The two models used for this study are wflow_sbm hydrologic model [22], and the existing Barotse hydrodynamic model [16] coupled using GLOFRIM v2.

Setting-up a Wflow_sbm v2019.2 Model for Upper Zambezi Basin
Wflow is a rainfall-runoff grid physically based distributed model, part of the Deltares open streams project. It uses open-source global data such as DEM, surface water network, land cover, soil types and their parameters to simulate hydrological processes from climate data such as precipitation, temperature, and radiation [22,23]. Its physical basis and use of global data make it suitable for cases where field observed data is lacking [22,24] and not adequate as is the case for the Upper Zambezi Catchment. To generate hydrologic input for this research, the hydrologic wflow_sbm at 2.15 arcmin spatial resolution (approximately 4×4 km at the Equator) was applied. The model was programmed in Python using the PCRaster Python extension [22]. The model consists of a set of python programs run on a command line to perform hydrological simulations over grid cells of static PCRaster maps [25,26]. Several water fluxes are calculated in wflow ( Figure 2). The Upper Zambezi Model has been set to run at daily time-steps from 2000 to 2018.
In this work, the model was forced with the fifth ECMWF ReAnalysis (ERA5) [27] for precipitation, potential evapotranspiration data and temperature and evaporation was computed using the Penman-Monteith equation. Globcover for landcover [28], SoilGrids [29] for soil parameters, we used MERIT Hydro; MERIT DEM for catchment delineation as well as Lake Hydro for lake data acquisition [30].

Setting Up Lisflood-FP v5.9.6 Model for a Subdomain Covering Barotse Sub-Catchment
In the simulation of wetland hydrodynamics, the existing Barotse Lisflood-FP Hydraulic Model was used [16]. The detailed characteristics and setup of Barotse Lisflood-FP Hydraulic Model used are not reported here but can be found in Makungu (2020) [16]. However, the Barotse Lisflood-FP Hydraulic Model was resampled to 500×500 meter in UTM 34S projection and was forced with runoff-on-a-grid as well as upstream inflows from WFLOW. Lisflood-FP is a two-dimensional (2D) programmed in Python C++, non-commercial raster-based flood inundation model designed for research purposes for simulating river flooding and floodplain inundation in data-scarce catchments by the University of Bristol. The hydraulic model solves numerically the local inertia [31], diffusive [32] or kinematic [33] approximations to the one-dimensional Saint-Venant equations to simulate the propagation of the flood wave through the river channel. In this model, the river channel is represented using the local inertia approximation implemented at sub-grid scale [31].

Setting Up a Barotse Floodplain Fully Coupled Hydrological-hydrodynamic Model Approach
We setup a fully coupled hydrological-hydrodynamic model framework using (wflow_sbm and Lisflood-FP), for the Barotse floodplains by adapting Hoch (2017a) framework ( Figure 3) [5]. This was done by defining the coupling settings in the GLOFRIM initialization file. This initialization file is then read by GLOFRIM to first initialize the model-specific configuration files and then initialized both models as a coupled unit. The coupling has been done per grid basis. Thus, once the two models are initialized, a loop is entered, beginning at the start time and terminating at the end time.

Wflow Sensitivity Analysis
We conducted sensitivity analysis on the WFLOW parameters with the sensitivity indices of P-value as shown in Table 2. Sensitivity analysis is the process of determining the rate of change in model output with respect to changes in model inputs (parameters). It is a necessary process in identify key parameters and parameter precision required for calibration [34]. The most sensitive parameters identified for use in the calibration of WFLOW model are as shown in Table 2. In choosing few parameters for calibration, we are conscious of the equifinality concept [35], that different combinations of values of the same parameters may produce an identical output signal.

RMSE-observations Standard Deviation
Ratio (RSR) The decision to use the combination of these four (4) metrics was guided by the recommendation of Ritter and Muñoz-CarSpena (2013) [37] to compensate weaknesses of some of the indicators. As they argued that when a single indicator is used it may lead to incorrect verification of the model. Instead, a combination of graphical results, absolute value error statistics and normalized goodness-of-fit statistics is ideal. The results of the model performance are also graphically shown in Figure 4 for calibration and Figure 5 for validation. Scholars have stated that a hydrological Model with the following matrices KGE >= 0.5, Pbias <= +/-25%, NSE >= 0.65, RSE >= 0.50 is generally considered to be a well performing model [37,38]. It can also be noted that the peak values for both observed and simulated are in similar ranges as are the low flows. However, looking at pbias, the model is slightly overestimating the flows in some years for both calibration and validation period bearing in mind that we used fifth generation (ERA5) of European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyses of the global climate for forcing data, which is not necessarily a good reproducer of absolute quantities of rainfall in comparison to absolute rainfall estimates from the local rainfall gauges. Nevertheless, for the purpose of this paper, the ERA5 forcing data was sufficient. As the focus of this work was on flow behavior and processes (interactions of groundwater and surface water) that we can capture using a fully coupled modelling.
Despite that wflow model's ability to simulate the upstream hydrology upstream very well, the wflow model results of the downstream Barotse Floodplain gauge station (Senanga) did not quite match the observed monthly flows as indicated by evaluation statistics: nse = 0.132; kge = 0.509; pbias = 37.740; RSR = 0.9233; and the hydrography ( Figure 6). This is likely because the representation of both floodplain channel hydrodynamics (storage, bifurcation, lateral connections) and vertical hydrological processes (floodplain water infiltration into the soil column; evapotranspiration from soil and vegetation; and evaporation of open water) is necessary to correctly capture floodplain dynamics as more processes are captured. This necessitates the need for a coupled model to capture feedback floodplain processes to improve floodplain simulation. This is evidenced from studies that have shown improvement in floodplain model simulations whenever more processes are captured by model coupling. For example, by using GLOFRIM, Hoch et al. (2019) [6] coupled the global hydrologic model PCR-GLOBWB with the hydrodynamic models CaMa-Flood and Lisflood-FP. Results show that replacing the kinematic wave approximation of the hydrologic model with the local inertia equation of CaMa-Flood greatly enhances the accuracy of peak discharge simulations, as expressed by an increase in the Nash-Sutcliffe efficiency (NSE) from 0.48 to 0.71. This is possibly because groundwater discharge to groundwater-dependent floodplain systems is a critical element in assessing such systems. Therefore, hydrological modelling alone is insufficient to capture significant processes such as groundwater discharge to the floodplain hydrodynamics. Hence, an approach that captures river flows, wetland inundation, and groundwater would be required to improve the assessment of floodplains for sustainable floodplain management.  [39] used the LISFLOOD hydrological model to simulate floodplain inflows and the Lisflood-FP hydraulic model to generate inundation extents. All these approaches are all based on one-way offline hydrologicalhydrodynamic model coupling. Although these approaches are useful for simulating natural or anthropic influences over wetland inflows, these methods are based on one-directional coupling, and they do not allow for two-way feedback between floodplain hydrodynamics and vertical hydrology (i.e., soil infiltration, evapotranspiration, etc.) as runoff generation and channel/floodplain dynamics are treated as independent processes. Consequently, as Hoch et al. (2017b) [18] argued, the resulting lack of feedback interactions between hydrology and hydrodynamics, for instance, by employing groundwater re-infiltration on inundated floodplains, can hamper modelled inundation and discharge results where such interactions may potentially be important, as in this case in the Barotse floodplain. Therefore, a two-way online fully model coupling approach allows for two-directional feedback between floodplain hydrodynamics and vertical hydrology since runoff generation and channel/floodplain dynamics are treated as interdependent processes. This is ideal for groundwater dependent wetlands such as the Barotse floodplain in order to capture vertical hydrological processes (baseflow contribution) and allow for a two-directional surface watergroundwater exchange.

Conclusion
To our knowledge, this is the very first attempt of promoting an approach for fully coupled model for Barotse Floodplain. A fully coupled model makes feasible a range of real-world problems that models based on kinematic or inertia-free approximations cannot simulate accurately. The exchange of variables between hydrology and hydrodynamics was done on a grid-to-grid basis at the time-step. The approach will serve as a basis for developing a fully coupled hydrological-hydrodynamic model for the Barotse floodplains. Fully coupled models have the potential to be used in implementing alternative wetland management strategies in the areas of water resources allocation, environmental flow (eflows), flood control, land use and climate change impact assessments, as well as pollution control. A fully coupled model coupling for Barotse presents an important step towards closing the gap between hydrology and hydrodynamics and for improved assessment of groundwater dependent wetlands in sub-Saharan Africa to inform the sustainable management of floodplain wetlands.

Data Availability Statement
Our coupled Model environment and approach for Barotse Floodplain are available at https://github.com/hcwinsemius/barotse.

Acknowledgements
We are grateful to the Worldwide Fund for Nature (WWF) Zambia and IWRM Center coordinators for supporting this work. Thanks, are also due to TU-Delft Faculty of Civil Engineering and Geosciences, Water Resources Section and Deltares for providing technical assistance. Thanks to the Department of Water Resources Development and Water Resources Management Agency (WARMA) for data on the Upper Zambezi Basin.

Institutional Review Board Statement
Not Applicable.

Informed Consent Statement
Not Applicable.

Declaration of Competing Interest
The authors declare that there is no conflict of interests regarding the publication of this manuscript. In addition, the ethical issues, including plagiarism, informed consent, misconduct, data fabrication and/or falsification, double publication and/or submission, and redundancies have been completely observed by the authors.