An application of the automatic domain updating to the Tonle Sap Lake, Cambodia

: Simulating flow dynamics in large-scale lakes is often time-consuming. For river flood simulation, the automatic domain updating (ADU), which can effectively control the simulation domain only in and around the flooded areas, has recently been developed. It is easily implementable without any computational errors for river flood simulation; however, its applicability to lake flow simulation with pre‐ cipitation/evapotranspiration has not been investigated. This study examines the applicability of the ADU to large-scale lake flow simulation with the 2-dimensional local inertial equations (2D-LIE) taking the Tonle Sap Lake, Cambodia, as a study site. The 2D-LIE with the ADU demonstrated 2.1 times faster simulation with errors less than 5.5%. This efficiency was achieved owing to the wet/dry seasonal nature of the tropical lake and backflow from the mainstream of the Mekong River in the rainy sea‐ son, suggesting that the ADU is applicable to large-scale lake flow simulation.


INTRODUCTION
Nowadays, assessing hydrodynamics of a freshwater body using mathematical models plays a fundamental role in designing lake conservation plans. In particular, largescale (> 100 km 2 ) lakes found in semi-arid areas, like the Tonle Sap Lake in Cambodia, are exposed to high risk of flood and/or contamination due to poor wastewater treatment (Ung et al., 2019;Hunsberger et al., 2018) in common with major megacities in Southeast Asia (Luo et al., 2019). The shallow water equations (SWEs) (Castro-Orgaz and Hager, 2019) is one of the most widely-used mathematical models for various large-scale water flow simulations such as river flooding (Luo et al., 2018) and/or coastal flooding (Ikeuchi et al., 2015). To reduce the computational burden of the SWEs, extensive collections of simplifications/approximations have been explored (Hunter et al., 2007;Teng et al., 2017;Mignot et al., 2019). One of the most effective approximations is the local inertial equations Correspondence to: Tomohiro Tanaka, Graduate School of Global Environmental Studies, Kyoto University, Room 584, C1-4, Kyotodaigakukatsura,Kyoto, (LIEs) first proposed by Bates et al. (2010), which is the SWEs without the advective acceleration terms (Seyoum et al., 2012;Martins et al., 2017;Mateo et al., 2017). This simplified model can overcome the severe numerical stability issue of the conventional diffusive wave equations, keeping a simpler form of the momentum equation than the SWEs (de Almeida and Yamazaki et al., 2013).
Another approach includes domain control techniques to reduce unnecessary processes in simulation. Yamaguchi et al. (2007) proposed the dynamic domain defining method. This method divides a simulation domain into several blocks and only computes the blocks containing wet cells to improve computational efficiency. A far simpler and more straightforward alternative, referred to as automatic domain updating (ADU), has recently been proposed by Tanaka et al. (2019). The ADU controls the simulation area at a cell level, i.e. it efficiently detects wet cells and their neighboring ones that may become wet in the next time step (shown as the blue and yellow cells in Figure 1, respectively). These domain control techniques can be easily implemented, reducing the computational costs by omitting the calculation of governing equations (and even omitting the conventional simple dry check process, i.e. checking if water depth is positive or not (Medeiros and Hagen, 2013) in dry areas: the green cells in Figure 1.
Since the ADU controls the simulation domain at a cell level in and around flood areas by tracing their boundaries at each time step, it implicitly assumes that the expansion of the flood areas starts only from specific (and thus known) cells. This assumption is reasonable when simulating river flooding events where the flood water always comes from river channels whose locations are known a priori. This scheme, in contrast, may not be transferable to simulating pluvial flooding events (floods caused by intensive rainfall exceeding drainage capacity) because the rainfall would be supplied to cells far away from the flood area as an external driver, and the flooding there cannot be straightforwardly handled by the ADU. Both pluvial and fluvial floods could happen in the same simulation domain. One example is large-scale lake flow simulation in a semiarid region where most of the heavy rainfall occurs in the rainy season during which the lake area widely expands (Kummu et al., 2014). In such a case, the impacts of omitting rainfall supply on dry cells might be limited because a large part of the target area is already wet (and included in the simulation area) in rainy seasons, while almost no rainfall is supplied in the dry seasons during which the lake shrinks.
This study examines the applicability of the ADU to a recently developed 2D lake flow model based on the LIEs on the Tonle Sap Lake (TSL) in the Mekong River basin, Cambodia, as a case study (Tanaka et al., 2018b). Computational results by the model with the ADU are compared with the one without the ADU focusing on the flood area and the observed data of the water stage. The impacts of ignorance of rainfall on dry cells, which is a key technique potentially contributing to efficient numerical computation in the ADU, are investigated as well. Finally, computational efficiency of the ADU is discussed for both parallel/ serial simulation settings assuming the different availability of computational resources. Figure 1. Definitions of the key terms and variables used in the ADU. The blue cells show the wet cells having water depths larger than a prescribed threshold value (0.01 (m) in this study); the yellow cells show the boundary cells; the green cells correspond to the dry cells. The numerical simulation takes place in total N wet wet (blue) cells and N b boundary (yellow) cells. If a yellow cell gets wet at a certain time step, it is updated as a wet cell (blue), and its surrounding dry cells are updated as boundary (yellow) cells. The momentum and mass balance equations are solved only at blue and yellow cells registered to an array C; the model checks dryness only at cells registered to an array B. We can accelerate the computation through limiting a process to check if a cell is dry or not only in the yellow cells

Governing equations
This study uses the lake flow model based on the 2D-LIE as the governing equations, developed by Tanaka et al. (2018b). The 2D-LIE consists of the continuity equation and the momentum equations where t is the time; z is the bed elevation, p is the line discharge in the x-direction, q is the line discharge in the ydirection, h is the water depth, r is the source term due to rainfall, e is the sink term due to evapotranspiration, Q tr is the inflow from tributary rivers, and g is the gravitational acceleration coefficient, and n is the Manning's roughness coefficient. The equations can be solved subject to appropriate initial and boundary conditions. Usually, the initial state of water depth and fluxes in the x-and y-directions are set as explained later; the water stage is prescribed at outflow boundaries, while the line discharges are specified at inflow boundaries.

Discretization
The 2D-LIE is discretized with a semi-implicit finite difference scheme (Bates et al., 2010), which is based on a spatio-temporally staggered discretization and an efficient treatment of the friction slope terms in the momentum Equations (2)-(3). The whole domain of water flows is uniformly discretized into an orthogonal structured computational mesh having rectangular cells. The time increment for the temporal integration is denoted as Δt. The spatial increments in the x and y directions are denoted as Δx and Δy, respectively. The quantity S evaluated at the time t = kΔt at the location (x, y) = (iΔx, iΔy) is denoted as S i, j k , where i, j, and k are integers. The quantities with halfinteger sub-and/or super-scripts are defined in the same manner.
Assume that ℎ i, j k , p i + 1/2, j k + 1/2 , and q i, j + 1/2 k + 1/2 have already been computed at each cell. We then discretize the 2D-LIE to temporally update these quantities in a staggered manner. Firstly, the continuity Equation (1) is discretized as Secondly, the momentum Equations (2)-(3) are discretized as

APPLICATION OF THE AUTOMATIC DOMAIN UPDATING
Here, the non-negativity preserving numerical water depths ℎ i + 1/2, j k and ℎ i, j + 1/2 k are evaluated as and respectively (de Almeida and . Although the friction term (the last term on the left hand side of Equations (5) and (6)) is discretized in a semi-implicit manner in terms of fluxes p i + 1/2, j k + 1/2 and q i, j + 1/2 k + 1/2 , they are expressed only by variables at the previous time step; therefore, it does not require any iterative procedures that conventional implicit schemes require. As suggested in the discretized equations, the water depth or the line discharge at each cell is calculated using the information coming from this and neighbor cells. Thereby, the ADU algorithm is easily implemented as detailed below.

Structure of the ADU
The ADU is a computational technique recently developed by Tanaka et al. (2019) for efficient river flood simulation. Theoretically, this method is applicable to any type of numerical simulation dealing with wave propagation phenomena expressed by hyperbolic partial differential equations with finite difference schemes. The core concept of the algorithm in the ADU is tracing the flood area and its surrounding cells by updating the boundary of the flood area at each time step (see also Figure 1).
The implementation procedure of the ADU is explained as follows: (1) To prepare the two vectors C = (c 1 , ..., c KL ) and B = (b 1 , ..., b KL ), representing the list of cells for calculation and for boundary updating, respectively. Here, K and L are the total number of cells in the x and y directions, respectively (2) To register N wet wet cells at the initial time to C and their surrounding N b cells to C and B (in the same way from initial time), where "wet" is defined as the status in which water depth is larger than a criterion (0.01 m in this study). This minimum water depth is set to avoid numerical instability caused by a too small water depth in Equations (8) and (9). The length of B is N b and that of C is (N wet + N b ) (3) To solve Equations (2) and (3) at all the cells belonging to C (4) To check all the cells in B; if b i becomes wet, to exclude it from the boundary list B. Then, add its surrounding cells to C and B (i.e. simulation domain expands and the boundary is updated accordingly). If it and all its surrounding cells become dry, the cell is excluded from B, and the surrounding ones not included in B are added to C and B (i.e. simulation domain shrinks), otherwise to include it in B again (5) To update the total number of boundary cells N b and the total number of cells for simulation (N wet + N b ) (6) To repeat Steps 3 and 4 until the terminal time of computation As indicated in the above procedures, this algorithm does not require any specific libraries. It is therefore transferable to models in any computer language. This study used a C++ language.

Implementation of ADU on lake flow simulation
The implementation of ADU to river flooding is largely different from that to lake flow simulations due to rainfall and evapotranspiration (right hand side of Equation (1)), while the latter is usually solely by the river discharge from the upstream. Because rainfall and/or evapotranspiration could happen at any cells over the domain, one may have to simulate the whole domain in case of not using the ADU. However, when it comes to floodplains having heavy rainfall and major floods in the same season, such as a semiarid area having dry and wet seasons, there is still room for applying the ADU to improve computational efficiency, by allowing rainfall and evapotranspiration only in wet cells as demonstrated later in this paper. These effects will be quantitatively evaluated in the following application.

Study site and simulation setting
The TSL is the largest lake in Southeast Asia, and is one of the largest aquatic environments for ecosystem exposed to social and climate changes (Salmivaara et al., 2016;Chan et al., 2019;Uk et al., 2018). The map of the study area is shown in Figure 2. The total area of the figure is 46,350 km 2 , in which the area surrounded by major highways with higher elevation (20,156 km 2 ) was set as the total simulation domain (where elevation data is shown in Figure 2). The actual simulation domain changes during the simulation by the ADU algorithm as explained above. Eleven tributary rivers flow into the TSL and lake water outflows into the Mekong River through the Tonle Sap River. The lake flow dynamics are strongly affected by the seasonal trend of the flow regimes in the Mekong River. In the rainy season from June to November, backflow occurs from the Mekong River to the TSL due to high water level, while the flooded water recedes to the Mekong River in the dry season from December to May.
The target period covers from July 1998 to December 2003 including a recent devastating flood event in 2000. Rainfall and temperature data are provided from gauge points (see Figure 2). Evapotranspiration is estimated from monthly temperature using the Thornthwaite method (Singh, 2016). Tributary river discharge data from the previous simulation result (Tanaka et al., 2018b) by the distributed hydrological model Geomorphology-Based Hydrological Model (Yang et al., 2002) is provided as the lateral inflow. The 1-D hydraulic model MIKE11 was used to generate the boundary water stage at the Prek Kdam station (Figure 2) (Fujii et al., 2003). Assuming the mild change in the water depth with a small flow velocity at the initial time of the simulation, the fluxes are set to be zero at the initial time. The initial water depth was set such that water stage becomes 1 m at grids with elevation lower than 1 m (corresponding to the main water body of the lake in the end of the dry season). Different and more realistic initial conditions can be set if such a condition is available based on some observation data, which was not available in our case. The 2D-LIE was numerically discretized with a spatial resolution of 500 m with an orthogonal structured computational mesh and a Manning's roughness coefficient of 0.03 m -1/3 s, which was used in Tanaka et al. (2018b) and confirmed that it reasonably reproduced the observed water stage at the Kampong Luuong station and outflow discharge at the Prek Kdam station.

Accuracy
The simulated water levels at the Kampong Luong station with/without the ADU are compared in Figure 3 (a). . Rainfall and temperature are observed at the points indicated by the light-blue and purple circles, respectively. Upstream tributary discharges (the green boxes) and downstream boundary water stage at the Prek Kdam station (the blue box) are given as the upstream and downstream boundary conditions, respectively. The data at the Kampong Luong station (the red box) was used for the validation of the 2D-LIE with the ADU They are so close to each other that the result with the ADU (red) is not visible during most of the simulation period (except for the dry season in 2001 and 2003), showing that the result with the ADU is only slightly smaller (0.15% in the dry season and 0.87% in the rainy season) than that without the ADU. Figure 3 (b) shows the time series of flood area in both simulations with/without the ADU (orange) and only in the simulation without the ADU (blue). A larger difference in the dry season appeared in remote areas far from the main lake where the water is supplied by local rainfall though its area is much smaller than the total flood area.
The average daily rainfall over the whole simulation domain used to drive the 2D-LIE with/without the ADU is  Figure 4 (recall that the ADU requires the rainfall data only in the wet cells), showing a clear difference over the simulation period (41.2% on average). The monthly mean of the missing ratio γ of the daily rainfall given to the 2D-LIE with the ADU (γ ADU ) to the one without it (γ Total ) is defined as follows ( Figure 5): As shown in Figure 5, the missing ratio goes up to 70% in the later dry season because of a larger amount of rainfall missing on the dry cells during the dry season when the lake, i.e. the simulation domain of the ADU, shrinks. On the other hand, their difference decreases in the later rainy season (October to November) to around 30%. The flood area expands the most in this season; thereby, the simulation domain also expands and the missing ratio becomes smaller. As the lake flow more dynamically changes during the rainy season, the above mechanism is considered to Figure 4. The time series of daily rainfall employed to externally drive the model with (red) and without (blue) the ADU. The ADU (red) missed some part of total rainfall given in the simulation without the ADU (blue) Figure 5. The monthly mean of the ratio of difference in daily rainfall between the 2D-LIE with the ADU to the one without the ADU. Missing ratio of rainfall in ADU is around 0.4 to 0.7 in the dry season and 0.1 to 0.3 in the rainy season contribute to the smaller impact of missing rainfall on simulation results. Furthermore, this impact on the water stage and/or flood area is particularly small in the TSL, owing to the dominant backflow from the Mekong River during the rainy season.

Computational cost
The 2D-LIE with/without the ADU was compiled with Intel@ C++ compiler 14 and carried out on a Linux machine with an Intel@ Xeon@ processor (CPU E5-2680 2.8 GHz). The computational time for the 5.5-year simulation from July 1998 to December 2003 was 43.0 and 91.7 hours with and without the ADU on a single core, respectively. The ADU thus accelerates the computation by a factor of 2.1. Tanaka et al. (2019) demonstrated 5 to 20 times faster simulation in a case study of river flooding in Japan, meaning that the computational efficiency of the ADU for the lake flow is less than those for river flooding. This large difference is simply explained by the ratio of flood area to the whole simulation domain. The flood area in simulating river flooding, particularly in the river confluence, has a variety of spatial distribution under difference scenarios. The flood area in such a case is thus significantly affected by the spatial distribution of rainfall and river defenses. On the other hand, the problem of the lake flow seems not to encounter this issue. Therefore, the simulation domain can be detected faster using the simple dry check. The improvement of the efficiency by the ADU then becomes much smaller. Computational time of each day during the simulation period is plotted with daily flood area in Figure 6. Their daily change is consistent with each other simply because flood area corresponds to the number of cells for computation. The effect of the ADU is clear in the dry season when flood area is smaller. The 2D-LIE with the ADU is 1.8 times and 1.4 times faster than without the ADU on 8 and 48 cores, respectively. Therefore, the ADU more significantly improves computational efficiency of the 2D-LIE as the computational environment becomes less massive. Our computational results suggest that the ADU is still an effective practical method for flood simulation involving Figure 6. Time series of computation time at each time step in the 2D-LIE with (red) and without (blue) the ADU and that of flood area (green). Computation time of each time step is consistent with flood area (green) for both simulations while its time is shorter in computation with the ADU (red) lake flow dynamics especially under limited computational resources. In summary, we find that the ADU algorithm enhances the numerical modelling of lake flow dynamics by significantly reducing the computational burden in exchange for the small reduction in accuracy in case of the large-scale lake flow simulation in semi-arid areas compared with conventional models.

CONCLUSIONS
The 2D-LIE enhanced with the ADU was applied to simulating surface water flows in and around the TSL in Cambodia. A theoretical estimate of the improvement of computational efficiency by the ADU was derived. The obtained computational results demonstrated that the total amount of rainfall used to drive the flow dynamics becomes much smaller with the ADU, mainly in the latter dry season and early rainy season during which the amount of rainfall is relatively small. Furthermore, the flow dynamics of the TSL in the rainy season is strongly affected by backflow from the Mekong River, which is qualitatively the same as river flooding. As a result, despite the significant ignorance of the rainfall in the model with the ADU, the 2D-LIE with/without the ADU are in good agreement in terms of the water stage and flood area. Note that this result is strongly affected by unique hydroclimatic features in large-scale lakes in semi-arid areas, and applicability of the ADU to areas driven by different mechanisms should be examined in future.
Future applications are expected in the assessment of the lake flow dynamics under climate change and/or anthropological impacts (e.g. dam operation in the Mekong River) scenarios, which require a computationally efficient model to complete a vast number of ensemble simulations (Liu et al., 2018;Tanaka et al., 2018a). The 2D-LIE with the ADU developed in this study can be used to address this topic. As the ADU itself is a technique to dynamically control the simulation domain, it is potentially applicable to various numerical models of flow simulations of importance in civil and environmental engineering. Such examples include the full shallow water equations, their non-Newtonian extension to debris flows (Ruiz-Villanueva et al., 2019;Han et al., 2019) and/or coupled solute transport phenomena (Liu, 2019), and multi-scale modelling (Li and Hodges;2019, Fukui et al., 2019. Future research would examine the ADU against these flow and transport phenomena to evaluate its computational performance in more detail.