The accurate and prompt mapping of flood-affected regions is important for effective disaster management, including damage assessment and relief efforts. While high-resolution optical imagery from satellites during disasters presents an opportunity for automated flood inundation mapping, existing segmentation models face challenges due to noises like cloud cover and tree canopies. Thanks to the digital elevation model (DEM) data readily available from sources such as United States Geological Survey (USGS), terrain guidance was utilized by recent graphical models such as hidden Markov trees (HMTs) to improve segmentation quality. Unfortunately, these methods either can only handle a small area where water levels at different locations are assumed to be consistent, or require restricted assumptions such as there is only one river channel. This paper presents an algorithm for flood extent mapping on large-scale Earth imagery, applicable to a large geographic area with multiple river channels. Since water level can vary a lot from upstream to downstream, we propose to detect river pixels in order to partition the remaining pixels into localized zones, each with a unique water level. In each zone, water at all locations flow to the same river entry point. Pixels in each zone are organized by an HMT to capture water flow directions guided by elevations. Moreover, a novel regularization scheme is designed to enforce inter-zone consistency by penalizing pixel-pairs of adjacent zones that violate terrain guidance. Efficient parallelization is made possible by coloring the zone adjacency graph to identify zones and zone-pairs that have no dependency and hence can be processed in parallel, and incremental one-pass terrain-guided scanning is conducted wherever applicable to reuse computations. Experiments demonstrate that our solution is more accurate than existing solutions, and can efficiently and accurately map out flooding pixels in a giant area of size 24,805 × 40,129. Despite the imbalanced workloads caused by a few large zonal HMTs dominating the serial computing time, our parallelization approach is effective and manages to achieve up to 14.3 × speedup on a machine with Intel Xeon Gold 6126 CPU @ 2.60GHz (24 cores, 48 threads) using 32 threads.