top of page

Chang Liao’s Research Showcase: Revisiting Flow Directions

I’ve been intrigued by modeling water flow since I earned my PhD in Earth and atmospheric studies. Much of my dissertation work focused on developing a three-dimensional model, ECO3D, of land-based ecosystems that I used to examine how lateral water flow impacts Arctic regions.

During my PhD defense someone asked me, “Can you apply your ECO3D model on the global scale?” (Liao, et al. 2019) This question inspired me to start thinking globally.

An extensive search of the literature showed no groups had performed watershed delineation, the basis of ECO3D, on a global scale. This led to work exploring how I could create a functional global model, which has taken me to some interesting intellectual places.

The main challenge to the above question is that Earth is not flat, meaning we cannot feed two-dimensional sections of the Earth into a traditional watershed delineation model.

So it’s reasonable to wonder how Earth scientists currently model global river flows. To date, over a dozen Global Hydrology Models (GHM) exist worldwide and all use latitude-longitude mesh grids to define global river networks. However, the latitude-longitude mesh grids in the Geographic Information System (GIS) are not all the same size (Figure 1).

Figure 1. Illustration showing the spatial distortion in latitude-longitude grids. Although the two sections are the same size according to the latitude-longitude grid, their actual geographic areas are substantially different.

This causes problems in defining flow directions, especially at high latitudes, as using grids with different sized squares introduces significant uncertainty to the surface runoff, a key hydrologic process. Can a different, non-latitude-longitude grid mesh create identically sized segments?

It can, and the answer lies in a familiar pattern—the soccer ball pattern of a hexagonal mesh, or Discrete Global Grid System (DGGS), can cover the Earth and facilitate uniform resolution (Figure 2).

Figure 2. Illustration of a Discrete Global Grid System that has uniform resolution (Sahr, 2015).

Shortly after joining Pacific Northwest National Laboratory (PNNL), I started working with other Earth scientists also interested in land surface and hydrologic modeling. Our group (Teklu Tesfa, Zhuoran Duan and L. Ruby Leung) started solving this problem from basic components: determining the appropriate tool for generating the hexagonal mesh, choosing the right programming language, selecting a testbed study area, etc.

It turns out that using a hexagonal mesh not only solves the problems associated with latitude-longitude grids but also provides other advantages and resolves several longstanding issues in computational hydrology. For example, it solves the so-called “diagonal flow path” and “island effect” issues (Johnston et al., 2009).

Designing a model implemented on a hexagonal mesh is just the beginning, though. As the hexagonal mesh differs dramatically from the traditional square grid structure, we had to rewrite all the algorithms used in the classical watershed delineation model. This process allowed me to develop a deep understanding of the fundamentals involved in terrain analysis. For example, I implemented the popular priority-flood algorithm that removes the local depressions on top of the hexagonal mesh (Figure 3) (Barnes et al., 2014).

Figure 3. Illustration of a priority-flood depression filling simulation using the hexagonal mesh method. Priority-flood is an efficient algorithm to fill DEM depressions by sequentially “flooding” the domain from a boundary inward to adjust elevations to assure surface drainage (Barnes et al., 2014). Light blue represents the initial default state; red represents a boundary; green represents a grid soon to-be-removed from the queue; orange represents a grid to-be-added into the queue; and purple represents finished grids. The numbers within a grid represent its global ID and elevation (in parentheses, unit: m), respectively.

To convert our work into something for wide use, we started at watershed scale and developed the freely available HexWatershed ( model. After evaluating and comparing the model to existing options, we found that using the hexagonal mesh-based flow routing could, if broadly applied, change the way the field approaches surface hydrology. For example, the new method not only provides the expected uniform flow direction (Figure 4), but also changes flow accumulation (Figure 5). This implies that the spatial pattern of evapotranspiration will change, meaning that all existing surface hydrology simulations could be improved.

Figure 4. The spatial distribution of simulated flow direction in a flat watershed in the Columbia River basin, WA, USA.

Figure 5. Histograms of flow accumulation from ArcSWAT (SGSD) and HexWatershed (HGSD) models, respectively. Between flow accumulation values of 0 and 4, the SGSD model consistently produces higher frequencies. In contrast, between flow accumulation values if 5 and 10, the HGSD model consistently shows higher frequencies.

Our team further applied the model at a global scale, attempting to tackle the original question that inspired my work. We produced the first uniform resolution global scale flow routing map, which can be used in a range of GHMs.

As often happens, answering one question leads to more questions. When I started the project, my primary interest was in terrestrial ecosystems. Meanwhile, PNNL’s Atmospheric Sciences & Global Change Division has scientists focused on the interactions between different Earth components (e.g., land-atmosphere, land-ocean, etc.). I soon realized the potential my work could have in coupled land-ocean simulations. HexWatershed is uniquely suited to couple land surface and oceanic models because the latter are usually based on unstructured meshes, like hexagonal grids.

This positions a further improved HexWatershed model to contribute to the United States Department of Energy’s Integrated Coastal Modeling (ICoM) project. One of ICoM’s key research focus areas uses a unified mesh, similar to a hexagonal mesh, to model the coastal zone continuum. The unified mesh framework aims to achieve the first seamless simulation of land-river-ocean interactions in areas affecting millions of people who live near coastal areas (Figure 6).

Figure 6. Illustration of current and next generation Earth system model (ESM) grid structures for the coastal zone continuum (Ward et al. 2020).

Developing HexWatershed is merely the first step in a long and continuing journey to create models that better represent water flows in a wide range of ecosystems.

I still plan on returning to my roots by applying the hexagonal framework to ECO3D or other GHMs to consistently model global terrestrial ecosystem dynamics, especially those in the climatically sensitive Arctic.


This research was partly funded by a Laboratory Directed Research and Development (LDRD) Program project and Earth System Model Development and Regional and Global Modeling and Analysis program areas of the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research as part of the multi-program, collaborative Integrated Coastal Modeling (ICoM) project.


Liao, C., Zhuang, Q., Leung, L. R., & Guo, L. (2019). Quantifying dissolved organic carbon dynamics using a three‐dimensional terrestrial ecosystem model at high spatial‐temporal resolutions. Journal of Advances in Modeling Earth Systems, 11, 4489– 4512.

Liao, C., Tesfa, T., Duan, Z., & Leung, L. R. (2020). Watershed delineation on a hexagonal mesh grid. Environmental Modelling & Software, 104702.

Barnes, R., Lehman, C., & Mulla, D. (2014). Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers & Geosciences, 62, 117–127.

Johnston, C. M., Dewald, T. G., Bondelid, T. R., Worstell, B. B., McKay, L. D., Rea, A., … Goodall, J. L. (2009). Evaluation of catchment delineation methods for the medium-resolution national hydrography dataset.

Sahr, K. (2015). DGGRID version 6.2 b: User documentation for discrete global grid software.

Ward, N.D., Megonigal, J.P., Bond-Lamberty, B. et al. Representing the function and sensitivity of coastal interfaces in Earth system models. Nat Commun 11, 2458 (2020).

By: Dr. Chang Liao

Earth Scientist

Pacific Northwest National Laboratory (PNNL)

7 views0 comments

Recent Posts

See All


bottom of page