Metadata-Version: 2.4
Name: baichuan-hydro
Version: 0.1.7
Summary: 百川（Baichuan）级联水文模型：多河道拓扑级联产汇流
Author-email: WLY <451215954@qq.com>
License: Proprietary
Project-URL: Homepage, https://example.com/baichuan
Project-URL: Documentation, https://example.com/baichuan
Project-URL: Source, https://example.com/baichuan
Project-URL: Issues, https://example.com/baichuan/issues
Keywords: hydrology,hydrological-model,runoff,routing,cascade,water-balance,水文,产汇流
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: Other/Proprietary License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: matplotlib>=3.5
Dynamic: license-file

# 百川（Baichuan）级联水文模型

> 百川：百川汇流、百川归海——多条河流（河道）逐级汇流成大河，恰如其分地描述本模型的
> **多河道拓扑级联汇流**核心。

逐小时、原创独立设计的**河道拓扑级联**水文模型。  
读取 OMS 格式参数文件（`params-testnew.csv`）+ 降雨文件（`data-test.csv`），模拟多河道的
**产流 → 河道内聚合 → 河道间拓扑汇流** 全过程，逐小时做水量平衡校验，并输出 `outlet.csv` / `QJ.csv`
（与 pyskby / Java 版布局一致）。通过 `BaichuanModel` 三步式 API 调用（见「快速开始」）。

---

## 版权与授权声明

**百川（Baichuan）** 为作者**从头独立设计**的简化水文模型，核心思想为「多河道拓扑级联汇流」。
本模型采用的均为**公开的水文学经典公式**（如 Hamon 蒸散公式、线性水库调蓄、水量平衡原理），
代码实现、数据结构、算法流程均为**独立原创**。

- 作者：WLY
- Copyright (c) 2026 WLY
- License：MIT（见 `LICENSE` 文件）
- 如有引用，建议标注本模型来源：「百川（Baichuan）级联水文模型」。

---

## 模型结构

## 模型结构

- **河道（HRU）= `@T nhru` 下每条记录**。每个河道有**独立出口流量**（即多个流量输出）。
- **河道拓扑由 `@T nchan` 段定义**：每条河道记录通过 `upst_inflow1` / `upst_inflow2` 列
  指明它接收哪些上游河道的来水（上游河道 `chan_id`，0 = 无）。例如 HRU2 的
  `upst_inflow1=1` 表示其上游是 HRU1 → 形成 HRU1 → HRU2 的级联。
- 每个河道（HRU）采用**等效单土壤模式**：`@T nhru` 直接提供 `field_capacity / conductivity /
  initial_water_content / porosity` 四个聚合列（等效土壤参数），不再逐土壤混合。
- 河道级调蓄参数 `chan_route_time`（滞留时间 τ，留存系数 `rk = exp(-1/τ)`）放在
  `@T nchan` 段（而非全局），每条河道可独立设置。

## 模型蓄量结构（四大水箱）

模型的水量平衡公式 `P = ET + Q + ΔS + ε` 中，蓄量变化 `ΔS` 由 **4 个水箱**的净变化构成：

| # | 水箱 | 蓄量变量 | 单位 | 进出水逻辑 |
|---|------|----------|------|-----------|
| 1 | **土壤水箱** | `hru["initial_water_content"]` | mm | 下渗水先入土；蓄满（>田间持水量 `field_capacity`）后多余部分溢出 |
| 2 | **壤中流水箱（SSR）** | `hru["ssstor_init"]` | mm | 土壤饱和溢出后暂存壤中流，按 `ssr2gw_rate` 缓慢退水 |
| 3 | **地下水水箱（GW）** | `hru["gwstor_init"]` | mm | 接受土壤补给 + SSR 补水，按 `gwsink_coef` 退水成基流 |
| 4 | **河道水箱（河道调蓄，多算法）** | `S_by_chan[cid]` | mm | 河道汇流演算，由 `chan_type` 选方法：默认线性水库 `outflow=(1-rk)(S+total_in)`（按 `chan_route_time` 留存），亦可选马斯京根-康吉(101)/运动波(1/3/4)/Lag(5) |

水量平衡校验项：`ΔS = ds(土壤) + dssr(壤中流) + dg(地下水) + dS_routing(河道)`。

> **水箱数量**：每个 HRU 有 3 个水箱（土壤、SSR、地下水），每个河道有 1 个路由水箱。
> 在"每河道 = 1 HRU"的模型里，2 河道示例共 `3×2 + 2 = 8` 个水箱。

### 四大水箱结构示意图

```
                    降雨 (precip)
                         │
                         ▼
   ┌─────────────────────────────────────────────┐
   │  ① 土壤水箱  initial_water_content (HRU级)  │
   │  下渗 → 蓄水；超过田间持水量 fc 后溢出      │
   │  蒸散 ET 从土壤水扣减                      │
   └──────────────┬──────────────────────────────┘
                  │ 超额水 excess（土壤蓄满后）
                  │
        ┌─────────┴──────────┐
        │                     │
        ▼                     ▼
   (按 soil2gw_max          ┌──────────────────────┐
     速率上限给地下)         │ ② 壤中流水箱 SSR    │
        │                   │ ssstor_init (HRU级)    │
        │                   │ 按 ssr2gw_rate 退水 │
        │                   └──────────┬───────────┘
        │                              │
        ▼                              ▼ 退出水按 ssr2gw 部分给地下
   ┌──────────────────────┐            │
   │ ③ 地下水水箱 GW      │◄───────────┘
   │ gwstor_init (HRU级)      │
   │ 按 gwsink_coef 退水  │──► 基流 baseflow ──┐
   └──────────────────────┘                    │
                                                ▼
   ┌──────────────────────────────────────────────┐
   │ ④ 河道水箱（线性水库） S_by_chan[cid]        │
   │ inflow = 本地产流 + 上游来水 + 外接水源       │
   │ outflow = (1-rk) × (S + inflow)              │
   └──────────────────┬───────────────────────────┘
                      │ 出口出流 outflow
                      ▼
              传给下游河道 / 汇出流域
```

#### 彩色结构示意图（PNG）

![四大水箱结构图](baichuan/structure.png)

> 该 PNG 由 `baichuan/plot_structure.py` 脚本生成（`python -m baichuan.plot_structure` 或 `baichuan/baichuan/plot_structure.py`）。

### 水量流动全链路

```
降雨 → 土壤水箱(①) →(蓄满溢出) 壤中流水箱(②) → 地下水水箱(③) → 基流
                    ↘ 土壤水直接 ET 蒸散
     壤中流出流(②)、基流(③)、超渗地表径流 → 河道水箱(④) → 出口出流 → 下游
```

## 参数文件结构（PRMS 风格）

`params-testnew.csv`（OMS `@T` 格式）含以下段落：

| 段 | 说明 |
|----|------|
| `@T nhru` | 每个 HRU = 一个河道：`hru_area, hru_psta, layer_depth, field_capacity, conductivity, initial_water_content, porosity, soil2gw_max`。`hru_psta` 为关联雨量站序号（从 1 开始，对应 `data-test.csv` 的 `precip[psta-1]` 列）；`layer_depth` 单位为 **m**（加载时 ×1000 转 mm）；`field_capacity`/`initial_water_content` 为体积含水率小数（× layer_depth 转蓄量 mm）；`conductivity` 为饱和水力传导度（**m/s**，加载时 ×3600×1000 转 mm/h） |
| `@T nchan` | **河道段**：`chan_id, chan_name, chan_route_time, upst_inflow1, upst_inflow2(3...)`。`chan_route_time` 为河道滞留时间 τ(h)，留存系数 `rk = exp(-1/τ)`；`upst_inflowN` 为上游河道引用（索引>0 才有效） |
| `@T ngw` | **地下水段**（PRMS 风格）：`ngw_id, gwflow_coef, gwsink_coef, gwstor_init, ngw_hru` |
| | - `gwsink_coef`：基流退水系数（= 旧 gw_recession），控制地下水出流快慢 |
| | - `gwstor_init`：地下水水箱初始蓄量（HRU 级） |
| | - `gwflow_coef`：保留字段（图片格式兼容），当前未参与计算 |
| | - `ngw_hru`：关联的 HRU id（一个 HRU 对应一个地下水单元） |
| | - 注：饱和超额的地下补给比例由 `soil2gw_max`（速率上限，mm/天）控制，不再用 `gw_split` 固定比例 |
| `@T nssr` | **壤中流（SSR）水箱段**（PRMS 风格）：`ssr_id, ssr2gw_rate, ssstor_init, ssr_hru` |
| | - `ssr2gw_rate`：SSR 退水系数（水箱放水快慢），同时是退出水给地下水的比例 |
| | - `ssstor_init`：SSR 水箱初始蓄量 |
| | - `ssr_hru`：关联的 HRU id（默认 1 对 1 对应 ngw，省略 ssr_gwres 字段） |

`data-test.csv`：OMS `@T obs` 逐小时降雨/观测文件，包含：

```
@T obs
date_start  yyyy MM dd HH mm ss
date_end    yyyy MM dd HH mm ss
date_form:yyyy MM dd HH mm ss
@H date,runoff[0],precip[0],precip[1],...
Type,date,real,real,real,...
2026 08 25 00 00 00,0.0,0.0,0.5
...
```

- 元数据行（`date_start` / `date_end` / `date_form`）仅说明时段与日期格式，程序读取时保留但不对时段做强制校验。
- `@H` 表头列：第一列为 `date`（时间戳），后续可含 `runoff[0]` 等观测径流列，以及 `precip[0]`, `precip[1]`, ... 等雨量站列。
- `Type` 行为类型说明行，程序自动忽略。
- 数据行第一列为时间戳，格式与 `date_form` 一致；后续列与 `@H` 对齐。
- 行数决定模拟时长。
- HRU 通过 `hru_psta` 关联雨量站：`hru_psta=1` 取 `precip[0]`，`hru_psta=2` 取 `precip[1]`，以此类推。

## 产汇流算法集成说明

百川将「**产流**」（HRU 级，降雨→径流）与「**汇流**」（河道级，河道间拓扑演算）解耦：
产流对所有河道统一采用等效单土箱 + 三层退水方案；汇流则按 `@T nchan` 的 `chan_type`
选择不同演算方法。**所有方法均满足逐小时水量平衡（连续方程）**。

### 一、产流方法（`calculate_hru_runoff`）

不区分 `chan_type`，每个 HRU 逐小时计算。采用**等效单土壤 + 三层线性水库退水**方案，
ET 从 HRU 总可用土壤水统一取（避免逐土壤重复扣减潜在 ET）。

| 子过程 | 算法名称 | 公式 | 物理意义 | 适用范围 |
|--------|---------|------|---------|---------|
| 潜在蒸散 PET | **Hamon 法** / **Epan 法**（二选一） | Hamon：`PET = f(日长, 气温)`；Epan：`PET = epan_depth × epan_coef` | 估算大气蒸发能力 | 缺实测水面蒸发用 Hamon；有蒸发皿观测用 Epan |
| 地表径流 | **超渗产流（Horton 型）** | `若 P > 下渗率 fc：Rs = P − fc；否则 Rs = 0` | 降雨强度超过饱和下渗率时产生 | 陡坡、低植被、暴雨初期 |
| 土壤填蓄 | **蓄满法（田间持水量）** | 下渗水优先填满 `fc − W`，填满后超额溢出 | 刻画土壤蓄水能力 | 所有下渗水 |
| 地下补给 | **速率上限法** | `gw_rech = min(soil2gw_max·Δt, excess)` | 限制每日补给速率，超出进壤中流 | 饱和后超额水拆分 |
| 壤中流 SSR | **线性水库（SSR 箱）** | `ssr_out = ssstor·ssr2gw_rate·(1−ssr2gw_rate)` | 壤中流产流 + 滞留 | 蓄满后侧向流 |
| 基流 | **线性退水（recession）** | `baseflow = gwstor·gwsink_coef` | 地下水缓慢退水成基流 | 枯水期基流 |

> **产流特点**：地表径流=超渗型（Horton）、壤中流+基流=蓄满型（Dunne），属「混合产流」；
> 全部输入来自参数文件，无额外计算量。若需纯蓄满产流（地表径流也来自饱和溢出）需改 `calculate_hru_runoff` 结构。

### 二、汇流方法（`chan_type` 分派）

`chan_type` 数字与 OMS/KINEROS2 保持一致；baichuan 自定义扩展从 **101** 起。
所有方法最终都等价于**拓扑级联 + 水量平衡**演算。

| chan_type | 算法名称 | 公式 / 递推 | 输入参数（读取） | 计算量（非输入） | 意义与适用范围 |
|-----------|---------|------------|----------------|----------------|--------------|
| **9**（默认） | **线性水库** | `rk = exp(−dt/res_k)`；`outflow = (1−rk)(S + total_in)` | `res_k` | 无 | 概念性调蓄，等价于 OMS "水库线性演算"；通用兜底 |
| **101** | **马斯京根-康吉 Muskingum-Cunge** | 先由几何推 K/X（`X = 0.5(1 − Q0·dt/(K·A0·c))`，波速 `c=1.67V`，`V` 由曼宁反算），再走马斯京根 `O = C0·I + C1·I_prev + C2·O_prev` | `chan_length`、`chan_slope`、`chan_rough`、`chan_width`、`chan_thresh`(起调流量 Q0) | **K、X**（Cunge 公式） | 物理性河道演算，无需率定 K/X；适用于有断面几何资料的天然/规则河道 |
| **1** | **矩形明渠运动波**（圣维南运动波，X=0 集总等价） | `K = L/c`，`c` 由 `Q=α·A^m` 关系得 | `chan_slope`、`chan_rough`、`chan_width`、`chan_thresh` | **α、m**（`ch_am` 矩形式：`α=1.49·√S0/(n·(width·3.28084)^0.6667)`，`m=1.67`）、`c`、`K` | 矩形断面河道；与 OMS/Java `ch_am` 严格数值对齐 |
| **3** | **三角明渠运动波**（同上，X=0） | 同 1 | `chan_slope`、`chan_rough`、`chan_t3_lbratio`、`chan_t3_rbratio`、`chan_thresh` | **α、m**（`ch_am` 三角式：`α=1.18·√S0/n·(√x1/(√x2+√x3))^(2/3)`，`m=1.33`）、`c`、`K` | 三角断面河道；与 OMS/Java `ch_am` 严格数值对齐 |
| **4** | **显式运动波参数** | `K = L/c`，`c = dQ/dA = α·m·A0^(m−1)` | `chan_alpha`、`chan_cmp`、`chan_thresh` | **c、K**（α/m 直接读，不反算） | 已知运动波关系 `Q=α·A^m` 的河道，最灵活 |
| **5** | **Lag 滞后模型**（纯平移） | `outflow[t] = inflow[t − Lag]`（线性插值） | `Lag`(=chan_lag) | 无（仅时间平移） | 忽略调蓄、只滞后；适用于短小陡急河道（**注：旁侧入流 rb_ssarea/rb_gwarea 暂未接入**） |
| **7** | **junction 汇流** | 仅上游组合（无调蓄） | 拓扑 `upstream` | — | 汇流节点；**当前未实现，走 9 兜底** |
| **8** | **水库 Puls 法** | 变库容曲线查表 `O = f(S)` | `res_low/res_up/res_mid` | — | 水库调蓄；**当前未实现（缺库容表），走 9 兜底** |
| **12** | **水库观测出流（RES_OBS）** | `outflow[t] = resObs[chan_resObs−1][t]`（缺失则 0）；`dS = total_in − outflow` | `chan_resObs`（指向降雨文件 `resObs[n-1]` 列）、降雨文件 `resObs[*]` 序列 | 无（出流为外部指定） | 标记为人造水库的河道：**跳过常规产汇流演算**，直接采用已知下泄过程，最贴近现实水库调度；上游来水仍透传下游 |

> **运动波（1/3/4）说明**：最终都集总等价成 **X=0 的马斯京根**（`route_muskingum(K, 0, …)`），
> 无扩散项；低于 `chan_thresh` 时 K→极小（近实时平移）。type=1/3 的 α/m 复刻 OMS/Java `ch_am`，
> 其中 **type=1 的 width 需 ×3.28084（米→英尺）** 以对齐 Java 英制公式常数 `1.49`。

### 二·五、人为干预：水库出流与分洪（叠加在常规汇流之上）

两类"外部指定"机制，均作用在**出流生成之后**，可独立或叠加使用，且与水量平衡完全兼容
（强制出流/分洪造成的净差额计入该河道蓄量 `dS_routing = total_in − outflow`，逐时闭合）：

| 机制 | 参数文件字段 | 降雨文件序列 | 公式 / 行为 | 适用范围 |
|------|------------|------------|-----------|---------|
| **水库观测出流** | `chan_resObs`（1-based） | `resObs[n-1]` | `type=12` 时 `outflow = resObs[chan_resObs−1]`，缺则 0；**跳过常规产汇流演算** | 仅 `chan_type=12` 河道 |
| **分洪控制** | `chan_fenhong`（1-based） | `fenhong[n-1]` | **所有 `chan_type` 生效**：`outflow = max(0, 常规出流 − fenhong[chan_fenhong−1])`，缺则 0 | 任意河道，含 type=12 |

> **分洪叠加顺序**：先由 `chan_type` 决定基础出流（含 type=12 的指定出流），**之后**再减去
> `fenhong` 序列。因此即使 type=12 已指定出流，仍可在其上再分洪，贴合"指定下泄 + 分洪调度"的现实。
> 被分走的水量同步从 `outlet.csv` 的 `fenhong[i]` 列输出（m³/s），供集成人员直接读取，无需回解析降雨文件。

### 三、输入参数分类

- **全部读取自参数文件 / 降雨文件**（无需计算）：产流全套 + 汇流各方法的几何/率定参数
  （`res_k`、`chan_length/slope/rough/width`、`chan_t3_lbratio/rbratio`、`chan_alpha/cmp`、
  `Lag` 等）；人为干预字段 `chan_resObs`、`chan_fenhong`（参数文件）及 `resObs[*]`、
  `fenhong[*]`（降雨文件观测/调度序列）。
- **计算量（非输入）**：汇流中的 `K`（蓄量常数）、`X`（权重，101 专用）、`α/m`（1/3 由 `ch_am` 算、
  4 直接读）、`c`（波速）、`rk`（退水系数，9 专用）、`PET`（Hamon/Epan）。

### 四、与 OMS/Java 的差异 & 未集成项

| 项目 | 状态 |
|------|------|
| 产流全套（等效土箱 + 三层退水） | ✅ 已实现 |
| 汇流 9 / 101 / 1 / 3 / 4 / 5 / 12 | ✅ 已实现；1/3 已与 Java `ch_am` 数值对齐；12 为水库观测出流（指定序列） |
| type=7 junction | ❌ 未实现，走 9 兜底 |
| type=8 水库 Puls 法 | ❌ 未实现（缺 `res_low/up/res_mid` 库容表），走 9 兜底 |
| type=5 旁侧入流（rb_ssarea/rb_gwarea） | ⚠️ 仅纯平移，未接 subsurface/groundwater 旁侧入流 |
| `chan_ndx`/`dx`/`dts`（有限差分分段） | ⚠️ 未用（百川为集总式，非有限差分） |

---

## 计算流程（A→B→C→D→E）

整个模拟是一个逐小时的时间步进循环。每时刻先按**拓扑排序**确定河道计算顺序（上游在前），
再对每个河道依次走完 A~E。

### A. 单 HRU 产流（`calculate_hru_runoff`，model.py:172）
对 HRU 内部的**每一种土壤**（占比 > 0）：
1. **蒸散 ET**：优先消耗该土壤水 `initial_water_content`；不够则蒸干（实际 ET = 剩余土壤水）。
2. **下渗与地表径流**：净雨 `net_rain` 超过下渗率 `conductivity` 的部分为地表径流，其余下渗。
3. **壤中流 / 地下补给**：下渗水若超过土壤蓄水能力（`field_capacity`），超额水按 `soil2gw_max`
   （速率上限，mm/天）限制进入地下补给 `gw_recharge`，其余进壤中流 `interflow`（守恒关键）。
4. **基流**：地下水为**河道（HRU）级**状态 `gwstor_init`（由 ngw 段 `gwstor_init` 初始化），
   所有土壤饱和后的地下补给汇入同一水库，按 `gwsink_coef`（= 旧 gw_recession）退水一次：
   `baseflow = gwstor_init × gwsink_coef`。

5. **壤中流水箱（SSR）**：土壤饱和后的**壤中流源**（= `excess − gw_recharge`）不再瞬时出流，
   而是先注入 `nssr` 段的 SSR 水箱（状态 `ssstor_init`，初值 `ssstor_init`）。每步水箱按
   `ssr2gw_rate` 退水：`ssr_release = ssstor_init × ssr2gw_rate`，其中
   `ssr2gw_rate` 比例转入地下水（`gwstor_init`），剩余 `1 - ssr2gw_rate` 作为壤中流出流进河道。
   这样既实现了壤中流的滞留（慢于地表径流、快于地下水），又让部分壤中水补给地下水。

→ 输出：单个河道（未经面积加权）的 `(q, et, dsoil, dgw, dssr)`，其中 `q` 含地表径流 + 壤中流 +
   河道级基流，`dgw`/`dssr` 分别为该河道地下蓄量、SSR 水箱的净变化（均参与水量平衡校验）。

#### 产流机制说明（地表径流 ≠ 土壤饱和溢出）

当前模型的**地表径流**与**壤中流/地下水补给**是**两件独立、分阶段**的事，不要混淆：

1. **地表径流（超渗产流 / Horton 型）**：在 ET 之后、下渗之前判定。
   净雨 `net_rain` 若超过下渗率 `conductivity`，超出部分直接成为地表径流，
   **此时水根本还没进土壤**。与土壤是否已经蓄满无关。
2. **壤中流 + 地下补给（蓄满产流 / Dunne 型）**：是**已经下渗进土壤的水**，
   在土壤蓄水 `initial_water_content` 达到 `field_capacity` 后，**装不下的 excess** 才按 `soil2gw_max`
   （速率上限）限制地下补给，其余进壤中流。

> 即：地表径流 = "雨下得太急、渗不下去"；壤中流/地下水 = "土已经喝饱、再也装不下"。

**举例**：某时刻 `net_rain = 10`，`conductivity = 5`，`initial_water_content` 距 `field_capacity` 还差 3，`soil2gw_max × td = 1`：
- 地表径流 = 10 − 5 = **5**（雨太急，一半没下渗，直接地面流走）
- 下渗 5：其中 3 存进 `initial_water_content`，剩余 excess = 2
- 地下补给 = min(soil2gw_max × td, excess) = min(1, 2) = **1**，壤中流 = 2 − 1 = **1**
- 最终该土壤出流 `q = 地表径流 5 + 壤中流 1 + 基流`

因此步骤③"土壤饱和后"那一步**只分流到壤中流与地下水**，地表径流在更早的超渗阶段已处理。
若需改为"纯蓄满产流"（地表径流也来自饱和后溢出），需调整 `calculate_hru_runoff` 结构。

> **等效单土壤模式**：`nhru` 提供 `field_capacity / conductivity / initial_water_content / porosity`
> 四列（等效土壤参数），模型把该 HRU 当作**一种等效土壤**直接计算，把加权好的等效值当作
> 单一土壤参数使用。

### B. 河道内聚合（`simulate` 内对单个河道调用 `calculate_hru_runoff`）
当前每个河道对应一个 HRU（`hru_area` 权重=1），故本地产流即 A 的结果；若未来一条河道下挂
多个 HRU 子单元，则在此按 `hru_area` 加权求和，合成该河道出口。同时把每个河道的逐时明细
存入 `hru_records` 供绘图。

- A 算的是"砖块"（子单元各算各的、零散结果）；
- B 把砖块砌成"这面墙（河道）"——得到每个河道的出口；
- 若一条河道只有一个子单元（权重=1），B 即平凡求和，看起来"A 直接给了出流"。

### C. 河道间拓扑汇流（`simulate` 内，按 `upst_inflow` 累加来水 + 河道调蓄）
对每个河道，用 `nchan` 段里的拓扑与参数做**河道与河道之间**的汇流：
- 本地产流 `q_local`（来自 B）；
- 叠加所有上游河道出口来水：`upstream_in = Σ outflow[up]`（up 取自 `upst_inflow1/upst_inflow2`）；
- 河道调蓄：由 `chan_type` 选择演算方法（详见上文「产汇流算法集成说明」），
  默认 `chan_type=9` 为线性水库（`chan_route_time` 为滞留时间 τ，`rk = exp(-1/τ)`：
  `outflow = (1 - rk) * (S + q_local + upstream_in)`）；亦可选 101 马斯京根-康吉、
  1/3/4 运动波、5 Lag 滞后等。

→ 这正是"河道间拓扑的汇流"所在，拓扑关系完全由 `@T nchan` 的 `upst_inflow` 列定义，
  不再依赖隐式的"上游=前半、下游=后半"拆分。

### D. 水量平衡校验（`simulate` 内，逐河道）
每小时对每个河道系统分别校验：
```
P + 上游来水 = ET + Q出口 + ΔS + ε
（P 本地降水，Q 出口出流，ΔS = Δsoil + Δssr + Δgw + ΔS_routing，ε 闭合误差）
```
输入 = 本地降水 + 上游来水；输出 = 出口出流 + ET；蓄变 = HRU 蓄变 + 河道蓄变。
`|ε| < 1e-6` 视为闭合（实测可达 1e-14 级浮点精度），否则打印 `[WARN]`。

### E. 输出与绘图（`export_csv` / `plot_hrus`）
- `export_csv`：逐小时结果写出 `model_output.csv`，按河道列出
  `precip_*/upstream_in_*/q_*/et_*/dsoil_*/dgw_*/drouting_*/outflow_*/bal_*`。
- `plot_hrus`：为**每个河道 HRU（HRU1/HRU2）** 各绘制一张独立过程线 PNG
  （产流 q / 蒸散 et / 土壤蓄变 dsoil / 地下蓄变 dgw），文件名为 `hru_HRU1.png` 等。

## 文件说明

| 文件 | 说明 |
|------|------|
| `baichuan/__init__.py` | 包入口：导出 `BaichuanModel`、`write_outlet_csv`、`write_qj_csv` 等 |
| `baichuan_model.py` | `BaichuanModel`：三步式封装（load → run → 写出） |
| `core.py` / `dataio.py` / `writers.py` | 计算核心 / 文件读写 / 结果写出（发布版为编译扩展） |
| `plot_hydrograph.py` | 产汇流过程线绘图模块（单图 / 所有河道子图） |
| `run.py` | 开发期运行入口（支持 groovy 案例，**不随 wheel 分发**） |
| `params-testnew.csv` | 参数文件（OMS `@T`：nhru / nchan / ngw / nssr 段） |
| `data-test.csv` | 逐小时降雨/观测文件（OMS `@T obs`，含 `precip[0]/precip[1]/...`，mm/h，行数=模拟时长） |
| `requirements.txt` | Python 依赖（matplotlib） |
| `LICENSE` | 授权协议（以 `pyproject.toml` 的 Proprietary 为准） |
| `README.md` | 本文档 |

## 快速开始

> **依赖**：`Python 3.9+`；绘图需 `matplotlib>=3.5`（仅绘图时需要，纯计算可不装）。

### 一、安装（发布版）

```bash
pip install baichuan-hydro==0.1.7
```

- 已构建 cp39-abi3 wheel（一份二进制在 CPython 3.9/3.10/3.11/3.12/3.13 通用）；其他 Python 版本从源码编译需本机有 Cython 与 C 编译器。
- **不要用 `0.1.0`**：该版本在 Windows 下传裸文件名（如 `"outlet.csv"`）会崩溃，已通过 `0.1.1` 修复。

### 二、Python API 三步式（推荐 · 第三方集成主接口）

模型对外只暴露一个对象 `BaichuanModel`，三步即可出结果：

```python
from baichuan import BaichuanModel

# params_path 必填；data_path 可省略（默认取参数文件同目录的 data-test.csv）
model = BaichuanModel(
    params_path="data/mian010705/params-testnew.csv",
    data_path="data/mian010705/data-test.csv",   # 可省略
    pet_method="epan",                           # "epan"(蒸发皿法,默认) 或 "hamon"(需气温)
)
model.run()
model.write_outlet_csv("outlet.csv", sim_name="mian010705")
model.write_qj_csv("QJ.csv",        sim_name="mian010705")
```

`run()` 之后可访问内部结果：`model.records`（逐时）、`model.topo_order`（拓扑序）、`model.chan_db`、`model.hrus_by_id`、`model.nchan`。
更低层函数也均可直接调用：`from baichuan import load_params, load_precip, simulate, write_outlet_csv, write_qj_csv`。

### 三、输出文件格式（OMS 兼容）

- **`outlet.csv`**：5 行头 + 数据行，字段含 `date, basin_ppt1, runoff[0], basin_cfs, LitleCon, Delay1..3, Fast1..3, basin_slowflow, basin_hortonian_mm, basin_ssflow_mm, basin_gwflow_mm, basin_actet, qcms[0..n-1], soil_moist_pct[0..n-1], qincms[0..n-1], hru_ppt[0..n-1], fenhong[0..n-1]`。
  - `qcms[i]`：河道 i **出流** (m³/s)；`qincms[i]`：河道 i **入流** (m³/s) = 上游来水之和 + 外水源；`basin_cfs`：流域出口流量 (m³/s)。
- **`QJ.csv`**：区间产流 `qcms_QJ[i]`，各河道**自身产流**(m³/s)，不含上游汇流。
- 两套输出与 pyskby / Java 版 `outlet.csv` / `QJ.csv` **列布局逐位一致**，下游成果展示工具无需改动即可直接解析。
- （另有 `model_output.csv` 为逐河道明细的**开发/诊断**导出，非对外交付格式。）

### 四、命令行（CLI）用法 —— 随 wheel 分发，pip 安装后即可用

`pyproject.toml` 的 `[project.scripts]` 已注册控制台入口 `baichuan = "baichuan.cli:main"`，
安装后终端直接键入 `baichuan`。它同时具备**计算**与**绘图**两类能力，底层都调用同一套
`BaichuanModel` 三步式 API（即上节的 Python 接口），不依赖开发仓的任何私有模块。

```bash
# 版本 / 帮助
baichuan --version
baichuan --help
baichuan run --help
baichuan plot --help
```

#### 4.1 运行模拟：`baichuan run`

两种调用方式，二选一：

**(A) groovy 兼容模式（最简，兼容旧 mian010705 调用习惯）**

只需给出 `sim` 名，CLI 会自动从 `<root>/sim/<sim>.groovy` 解析出 `parameter(file:)`
指向的参数文件与 `inputFile` 指向的降雨文件，跑完写出结果：

```bash
# --root 指向"项目根"，即同时含有 sim/ 与 data/ 的目录（默认=当前目录）
baichuan run mian010705 --root /path/to/project_root
```

- groovy 文件位置：`<root>/sim/mian010705.groovy`
  （例如 mian010705 项目里它引用 `data/mian010705/params-testnew.csv` 与 `data/mian010705/data-test.csv`）
- 输出位置：`<root>/output/mian010705/out/outlet.csv` 与 `QJ.csv`
- 若 `--root` 下找不到 `sim/<sim>.groovy`，CLI 会打印带 `--root` 使用提示的
  `FileNotFoundError` 并退出码 2，**不会**误用默认/其他数据。

> **关于 `--root` 默认值的重要提醒**：省略 `--root` 时，它取**执行命令时 shell 所在的当前工作目录**（即你在哪个目录下敲的 `baichuan`，root 就是哪个目录），**既不是包所在目录，也不是 `site-packages/`**。
> 因此客户 `pip install` 后，包装在 `site-packages/`，但 `--root` 默认是**调用位置**的当前目录，跟包在哪无关。客户必须**自己提供**含 `sim/` 与 `data/` 的项目根——要么先 `cd` 进自己的项目根再执行 `baichuan run mian010705`，要么显式加 `--root /path/to/their/project`。
> 不要指望"装好包就自带 `sim/`"；`sim/` 与 `data/` 需随项目一并交付客户，或改用 `(B)` 显式 `--params/--data` 模式彻底不依赖 groovy。

**(B) 显式文件模式（不依赖 groovy，推荐给新接入的第三方）**

```bash
baichuan run \
  --params data/mian010705/params-testnew.csv \
  --data   data/mian010705/data-test.csv \
  --out    outlet.csv \
  --qj     QJ.csv
```

- `--params` 必填；`--data` 省略时默认取参数文件同目录下的 `data-test.csv`
- `--out` / `--qj` 省略时分别默认当前目录的 `outlet.csv` / `QJ.csv`

**通用可选参数**（A/B 模式都支持）：

| 参数 | 说明 |
|------|------|
| `--pet {epan,hamon}` | 潜在蒸散算法：`epan` 蒸发皿法（默认）/ `hamon` 需气温 |
| `--plot` | 计算完成后自动绘制所有河道产汇流过程线 PNG |
| `--out/--qj` | 自定义 outlet / QJ 输出路径 |
| `--sim-name` | 写入结果表头的 sim 名（默认从 sim 名或输出文件名推导） |
| `--root` | groovy 模式下的项目根目录（含 `sim/` 与 `data/`）；**省略则默认命令调用处的当前目录**（非包目录） |

#### 4.2 仅绘图：`baichuan plot`

复用已模拟结果或直接加载示例文件绘图（不跑模拟）：

```bash
baichuan plot --all            # 一个文件画出所有河道（每河道一个子图）
baichuan plot --chan 1         # 指定河道（hru_id，默认=出口河道）
baichuan plot --params P.csv --precip D.csv --out fig.png   # 自定义输入/输出
```

- 默认 `--params` / `--precip` 为包内示例 `model_params.csv` / `model_precip.csv`。

> 控制台运行时会逐时刻打印各河道出口流量，结束时给出总步数、出口流量与平均水量平衡
> 闭合误差（`1e-14` 量级，浮点精度内即视为守恒）。

### 五、开发期命令（不随 wheel 分发）

- **`python run.py <sim> --root <项目根>`**（开发示例脚本，**不随 wheel 分发**）：以 `sim/<sim>.groovy` 为入口，从 groovy 内 `parameter(file:)` / `inputFile` 解析参数与降雨 CSV，跑完写出 `output/<sim>/out/outlet.csv` 与 `QJ.csv`。
  ```bash
  python run.py mian010705                 # 默认 --root 为 baichuan/ 的父目录（项目根）
  python run.py mian010705 --pet-method hamon
  ```
  - 它与 `baichuan run <sim>` 在功能上等价，但 `run.py` 不在包内、强依赖开发仓的 `sim/` 与 `data/` 布局，**不应发布给客户**。面向客户的兼容层请用上节的 `baichuan run <sim>`。

### 六、输入文件约定

输入为 **OMS `@T` 格式**的参数文件 + 降雨文件（与 pyskby / Java 同源），以 `data/mian010705/params-testnew.csv` 与 `data/mian010705/data-test.csv` 为模板即可：

| 文件 | 格式 | 关键段落 / 列 |
|------|------|---------------|
| 参数 `params-testnew.csv` | `@T` 多段 | `@T nhru`（每 HRU=一河道）、`@T nchan`（拓扑 `upst_inflow*`）、`@T ngw`、`@T nssr` |
| 降雨 `data-test.csv` | `@T obs` | `date` + `precip[0],precip[1],...`（mm/h，行数=模拟时长）；HRU 经 `hru_psta` 关联 `precip[hru_psta-1]` |

段落/列定义详见下文「输入 / 输出参数清单」。新增河道只需在 `@T nhru` / `@T nchan` 各加一行并对应加 `precip[]` 列，模型自动拓扑排序，无需改代码。

### 七、注意事项

- 输出路径传**裸文件名**（如 `"outlet.csv"`）也能正常工作（`0.1.1` 已修复 Windows 下 `makedirs` 崩溃）。
- 逐时水量平衡闭合误差约 `1e-14`，量级为浮点精度，模型守恒。
- **授权**：`pyproject.toml` 标注为 `Proprietary`（以随包 `LICENSE` 为准），分发给第三方前请与作者确认授权范围。



### 八、打包发布流程（维护者）

发布到 PyPI 的标准步骤（以升级到新版本为例）：

1. **升版本号**：同步改四处
   - `pyproject.toml` 的 `version`
   - `baichuan/__init__.py` 的 `__version__`（CLI 的 `--version` 主读安装元数据，这里保持同步）
   - `baichuan/cli.py` 中 `--version` 的回退版本字符串（`except` 分支里的 `_pkg_ver`）
   - `README.md` 里 `pip install baichuan-hydro==X.Y.Z` 的版本
2. **构建 wheel**（只构建 wheel，**不要 sdist**，避免源码外泄）：
   ```bash
   python build_dist.py
   ```
   该脚本会先清 `build/`、`dist/`、`*.egg-info/` 再 `python -m build --wheel --no-isolation`。
   **清理这步不能省**：受保护模块的 `.py` 是在构建期由 `setup.py` 的 `build_py` 子类挡下的，
   而上一轮遗留在 `build/lib` 里的 `.py` 会被 `bdist_wheel` 原样打进 wheel。
3. **（推荐）本地验证**：装进临时 venv 跑 `baichuan --version` 与一个 `baichuan run`，确认 `.pyd` 可独立运行。
4. **上传 PyPI**（仅传 wheel；twine 进度条在本机 GBK 控制台会崩，需关闭并设 UTF-8）：
   ```bash
   $env:PYTHONIOENCODING="utf-8"
   python -m twine upload --disable-progress-bar dist/baichuan_hydro-X.Y.Z-*.whl
   ```

> **为什么这样排除源码**：`pyproject.toml` 的 `exclude-package-data` 只能排除*非 Python 的包数据*（可删 `*.c`），对 `.py` 模块无效；顶层 `exclude` 键在 setuptools 里也不合法。因此由 `setup.py` 的 `build_py` 子类在构建期过滤 `find_package_modules()`，被 Cython 编译的核心模块（`core`/`dataio`/`baichuan_model`/`model`/`writers`）的同名 `.py` 根本不进 `build/lib`，wheel 从生成那一刻起就只含 `.pyd`，无需事后拆包。

## 扩展：增加河道

1. 在参数文件（如 `params-testnew.csv`）的 `@T nhru` 增加一行（`hru_area, hru_psta, layer_depth, field_capacity, conductivity, initial_water_content, porosity, soil2gw_max`）。
2. 在 `@T nchan` 增加一行，设置 `chan_route_time`（滞留时间 τ）与 `upst_inflow1/upst_inflow2(3...)`（指向其上游河道 id）。
4. 在降雨文件（如 `data-test.csv`）的 `@H` 行增加一列 `precip[N]`（N 从 0 开始连续），并将该 HRU 的 `hru_psta` 设为 `N+1`。
5. 在数据行追加对应 `precip[N]` 的逐时降雨量。

模型会自动完成拓扑排序并级联计算，无需改代码。

---

## 输入 / 输出参数清单

### 一、输入参数

#### 1. 河道 / HRU 定义参数（`@T nhru`，每 HRU = 一个河道）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `hru_area` | 面积权重 | 该河道控制面积权重（加权汇流用），当前均=1 |
| `hru_psta` | — | 关联雨量站序号（从 1 开始，对应 `data-test.csv` 的 `precip[psta-1]` 列） |
| `layer_depth` | **m** | 土层深度（加载时 ×1000 转 mm），将体积含水率换算为蓄水量(mm) |
| `field_capacity` | 体积含水率 | 等效田间持水量（小数 0~1），× layer_depth(m→mm) → 蓄量 mm |
| `conductivity` | **m/s** | 饱和水力传导度/下渗率（加载时 ×3600×1000 转 mm/h），净雨超过部分为地表径流（超渗产流） |
| `initial_water_content` | 体积含水率 | 初始土壤含水率（小数 0~1），× layer_depth(m→mm) → 初始蓄量 mm |
| `porosity` | 体积比 | 孔隙率（小数） |
| `soil2gw_max` | mm/天 | 土壤→地下水补给速率上限 |

#### 2. 河道拓扑参数（`@T nchan`，每河道一行）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `chan_id` | — | 河道编号（与 `hru_id` 对应） |
| `chan_name` / `rvcd` | — | 河道名称（输出列名前缀） |
| `chan_route_time` | h | 河道滞留时间 τ，留存系数 `rk = exp(-1/τ)` |
| `upst_inflow1` / `upst_inflow2` / `upst_inflow3`... | — | 上游河道 `chan_id`（索引>0 才有效，0=无），有多少个就解析多少个 |
| `chan_source1` / `chan_source2` | — | 外接水源引用（`source[]` 数组的 1-based 索引，>0 有效） |

#### 3. 地下水参数（`@T ngw`，每 HRU 一行）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `ngw_hru` | — | 该地下水单元关联的 HRU id |
| `gwflow_coef` | 日退水系数 | 地下水→河道退水系数（日量级，由 `daily_to_hourly` 折算逐时） |
| `gwsink_coef` | 日退水系数 | 基流退水系数（= 旧 gw_recession） |
| `gwstor_init` | mm | 地下水水箱初始蓄量 |

#### 4. 壤中流（SSR）水箱参数（`@T nssr`，每 HRU 一行）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `ssr_hru` | — | 该壤中流水箱关联的 HRU id |
| `ssr2gw_rate` | 日退水系数 | SSR 退水系数（水箱放水快慢），同时是退出水给地下水的比例 |
| `ssstor_init` | mm | SSR 水箱初始蓄量 |

#### 5. 地下水参数（`params-testnew.csv` → `@T ngw`，每 HRU 一行）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `gwflow_coef` | 日退水系数 | 地下水→河道主退水系数（日量级，由 load_params 折算为逐时当量；当前主退水走 gwsink_coef，本字段预留） |
| `gwsink_coef` | 日退水系数 | 基流退水系数（=旧 `gw_recession`，日量级折算逐时），`baseflow = gwstor_init × gwsink_coef` |
| `gwstor_init` | mm | 地下水水箱初始蓄量（HRU 级） |
| `soil2gw_max` | mm/天 | 土壤→地下水补给速率上限（速率上限），控制饱和超额中进地下补给的比例：`gw_recharge = min(soil2gw_max × td, excess)` |

#### 6. 壤中流（SSR）水箱参数（`params-testnew.csv` → `@T nssr`，每 HRU 一行）

| 参数名 | 单位 | 含义 |
|--------|------|------|
| `ssr2gw_rate` | — | SSR 退水系数（水箱放水快慢），同时是退出水给地下水的比例 |
| `ssstor_init` | mm | SSR 水箱初始蓄量 |

#### 7. 降雨输入（`data-test.csv` → `@T obs`）

| 列名 | 单位 | 含义 |
|------|------|------|
| `date` | — | 时刻（格式由 `date_form` 指定） |
| `runoff[0]` | mm/h | 观测径流列（保留字段，当前未参与计算） |
| `precip[0], precip[1], ...` | mm/h | 各雨量站逐小时降雨量；`precip[N]` 对应 `hru_psta = N+1`，行数=模拟时长 |

---

### 二、输出参数

#### 1. 控制台打印

- 逐时刻各河道 **出口流量** `outflow`（mm/h，保留 3 位小数）
- 逐河道 **水量平衡闭合误差** `bal`（科学计数法，1e-6 内视为闭合）

#### 2. `model_output.csv`（按河道列出，N=HRU1/HRU2...）

| 列名 | 单位 | 含义 |
|------|------|------|
| `timestamp` | — | 时刻 |
| `precip_{N}` | mm/h | 该河道本地降水 |
| `upstream_in_{N}` | mm/h | 上游河道汇入来水（=Σ 上游 `outflow`） |
| `q_{N}` | mm/h | 本地产流（地表径流 + 壤中流 + 河道级基流，未经面积加权） |
| `et_{N}` | mm/h | 实际蒸散 |
| `dsoil_{N}` | mm | 土壤蓄量净变化 Δsoil（per-soil 加权） |
| `dgw_{N}` | mm | 河道级地下蓄量净变化 Δgw |
| `dssr_{N}` | mm | 壤中流（SSR）水箱蓄量净变化 |
| `drouting_{N}` | mm | 河道线性水库蓄量净变化 ΔS_routing |
| `outflow_{N}` | mm/h | 该河道最终出口流量（产汇流后） |
| `bal_{N}` | mm/h | 水量平衡闭合误差：`P + 上游来水 − (ET + Q出口 + Δsoil + Δgw + Δssr + ΔS_routing)` |

#### 3. 过程线图（`hru_HRU1.png` / `hru_HRU2.png` 等）

每张图自上而下 5~6 个子图（单位均为 mm/h 或 mm）：

| 子图 | 变量 | 含义 |
|------|------|------|
| ① | `outflow` | 出口流量（产汇流后最终出口） |
| ② | `q` | 本地产流 |
| ③ | `et` | 蒸散 |
| ④ | `dsoil` | 土壤蓄变 |
| ⑤ | `dgw` | 地下蓄变 |
| ⑥ | `dssr` | 壤中流水箱蓄变（有 SSR 时附加） |

> 注：所有"蓄变"类（dsoil/dgw/dssr/drouting）为**时段净变化量（mm）**，其余通量为 **mm/h**。

### 三、参数单位汇总

| 类别 | 单位 |
|------|------|
| 通量（降水/产流/蒸散/出口） | mm/h |
| 蓄量 / 蓄变（initial_water_content/field_capacity/gwstor_init/ssstor_init 及 Δ） | mm |
| 长度速率上限（soil2gw_max） | mm/天（速率上限，控制饱和超额中进地下补给的比例） |
| 退水系数（ssr2gw_rate/gwsink_coef/gwflow_coef） | 无量纲（0~1，日退水系数，由 `daily_to_hourly` 折算逐时） |
| 滞留时间（chan_route_time） | h（留存系数 rk = exp(-1/τ)） |
| 面积权重 / 土壤占比 | 无量纲（0~1，累加=1） |

> **单位说明**：本模型整体采用 **毫米(mm)**，参数文件统一按 **OMS 单位约定**：
>
> | 参数 | 参数文件单位 | 内部换算 |
> |------|-------------|----------|
> | `layer_depth` | **m** | ×1000 → mm |
> | `conductivity` | **m/s**（饱和水力传导度） | ×3600×1000 → mm/h |
> | `field_capacity` / `initial_water_content` | 体积含水率（小数 0~1） | × layer_depth(m→mm) → 蓄量 mm |
> | `soil2gw_max` | mm/天（土壤→地下水补给速率上限） | — |
> | `gwflow_coef` / `gwsink_coef` / `ssr2gw_rate` | 日退水系数（无量纲） | `daily_to_hourly` 折算逐时当量 |
> | `chan_route_time` | 小时 h（河道滞留时间 τ） | `rk = exp(-1/τ)` 留存系数 |
> | `hru_area` | 面积权重（无量纲） | — |
>
> 对应地下补给计算为 `gw_recharge = min(soil2gw_max × td, excess)`，其中 `td = deltim/24`（天）。
> 若与其他以 **英寸(in)** 为单位的模型对比，仅需长度换算（1 in = 25.4 mm）。

---

## 引用格式

如在论文或项目中引用本模型，建议采用以下格式：

> WLY. 百川（Baichuan）级联水文模型[CP/OL]. 2026. MIT License.

或英文：

> WLY. Baichuan: A Cascade Channel Hydrological Model[CP/OL]. 2026. MIT License.

## 文件清单

| 文件 | 说明 |
|------|------|
| `baichuan/__init__.py` | 包入口：导出 `BaichuanModel`、`write_outlet_csv`、`write_qj_csv` 等 |
| `baichuan_model.py` | `BaichuanModel`：三步式封装（load → run → 写出） |
| `core.py` / `dataio.py` / `writers.py` | 计算核心 / 文件读写 / 结果写出（发布版为编译扩展） |
| `plot_hydrograph.py` | 产汇流过程线绘图模块（单图 / 所有河道子图） |
| `run.py` | 开发期运行入口（支持 groovy 案例，**不随 wheel 分发**） |
| `params-testnew.csv` | 参数文件（OMS `@T`：nhru / nchan / ngw / nssr 段） |
| `data-test.csv` | 逐小时降雨/观测文件（OMS `@T obs`，含 `precip[0]/precip[1]/...`） |
| `requirements.txt` | Python 依赖（matplotlib） |
| `LICENSE` | 授权协议（以 `pyproject.toml` 的 Proprietary 为准） |
| `README.md` | 本文档 |
