Flopy/MODFLOW的DISU包是否包含单元格ID对应的层标识符？

架构老林

2026-7-18

在Flopy的DISU包中关联单元格ID与层号的方法

DISU包本身并没有内置存储单元格ID到层号的直接映射——因为非结构化网格（DISU）的设计本身不强制“层”的结构，单元格仅通过索引标识。但可以通过以下两种通用方法实现关联，满足你仅估算第二层渗透系数的需求：

建模阶段主动维护映射关系
在创建DISU包时，如果你是按层的顺序定义单元格（这是常见的建模方式），可以同步生成一个cell_to_layer数组，数组索引对应单元格ID，数组值对应所属层号。比如：
```
nlay = 3  # 模型总层数
ncells_per_layer = [100, 150, 120]  # 每层单元格数量
cell_to_layer = []
for layer_idx in range(nlay):
    # 给当前层的所有单元格标记对应的层号（比如从1开始计数）
    cell_to_layer.extend([layer_idx + 1] * ncells_per_layer[layer_idx])
```
之后要筛选第二层的单元格，只需执行：
```
layer_2_cell_ids = [cid for cid, lyr in enumerate(cell_to_layer) if lyr == 2]
```
这种方式完全通用，只要建模时保持单元格按层组织的逻辑即可。
通过单元格几何信息反向推导
如果没有提前维护映射，可以利用单元格的顶底高程信息分组。DISU的top和botm属性存储了每个单元格的顶、底高程，同一层的单元格通常具有相近的垂直范围。你可以通过统计高程区间，将单元格划分到对应层中，比如：
```
# 获取所有单元格的顶高程
cell_tops = disu.top.array
# 假设层的分界高程已知（比如layer_breaks = [100, 80, 60, 40]）
layer_breaks = [100, 80, 60, 40]
cell_to_layer = []
for top in cell_tops:
    # 判断当前单元格属于哪个层
    for lyr in range(1, nlay+1):
        if layer_breaks[lyr-1] > top >= layer_breaks[lyr]:
            cell_to_layer.append(lyr)
            break
```
这种方式适合已有的DISU模型，无需修改建模流程，但需要明确层的高程分界。