Chapter 3 SHUD模型软件与数据

3.1 安装 SHUD and rSHUD

本章介绍SHUD模型和rSHUD的安装过程,以及必要的开发库的安装。

3.1.1 SUNDIALS/CVODE

SHUD模型的求解需要SUNDIAL中CVODE库的支持。 SUNDIALS(SUite of Nonlinear and Differential/ALgebraic equation Solvers)由六个求解器组成;其中的 CVODE 是一个用于求解刚性和非刚性常微分方程系统(初始值问题)的数值求解器,数学表达为: \(y' = f(t,y)\). CVODE支持可变阶精度、可变步长数值方法。CVODE字母组合含义为“C版本的可变系数常微分方程求解器(Variable-Coefficient Ordinary Differential Equation Solver)” 。使用SHUD模型只需要CVODE支持,因此无需安装SUNDIALS所有求解器,只需要安装CVODE即可。

因为SUNDIALS/CVODE持续更新中,代码和函数定义等等变化较大,因此SHUD模型通常仅支持CVODE较新版本,并非所有版本都可支持。

SUNDIALS/CVODE的安装方法有两种:(1)自行下载安装包; (2)使用SHUD源代码的安装脚本。

3.1.1.1 使用脚本安装CVODE

源代码文件夹中的configure.sh文件是自动安装CVODE的脚本——可以自动从github上下载CVODE最新源代码并通过cmake编译并安装在用户目录下。

3.1.1.2 自行下载并安装CVODE

SUNDIALS/CVODE 代码的下载地址为:https://computation.llnl.gov/projects/sundials/sundials-software

CVODE 5.x安装步骤:

  1. 进入命令行界面,并解压下载好的cvode压缩包。
  2. 在CVODE的根目录下新建文件夹 builddir, 并进入该文件夹。
mkdir builddir
cd builddir/
  1. 测试 ccmake. 如果没有 cmake请安装。
ccmake 
  1. 运行ccmake配置安装环境。
ccmake ../sundials/cvode-5.0.0
Screenshot of Step 1 in ccmake
Screenshot of Step 1 in ccmake

首先出现的是空白配置页面。按c开始配置。

Screen short of Step 2 in ccmake
Screen short of Step 2 in ccmake

上图为默认的配置信息图片上的四处配置信息需要注意,值需要改为下面的值。OPENMP值需要您的系统支持OpenMP,因此在确认您的系统支持OpenMP的情况下打开。

BUILD_CVODE = ON
CMAKE_INSTALL_PREFIX = ~/sundials
EXAMPLES_INSTALL_PATH = ~/sundials/examples

以上值修改完成后,按c确认配置。

Screenshot of Step 3 in ccmake
Screenshot of Step 3 in ccmake

ccmake开始自动配置安装参数。当配置完成后,按g生成配置文件并退出。

  1. 然后使用make编译并安装:
make
make install 

3.1.2 编译SHUD

编译信息保存在文件Makefile,其中关键的编译参数:

  1. SUNDIALS_DIR的保存路径(非常关键). 默认值为 ~/sundials,如果您已经安装到此路径,即可保持默认值。
  2. OpenMP的路径。当您的系统支持OpenMP时,可根据实际安装路径配置。
  3. SRC_DIR路径。默认值为 SRC_DIR = .,即当前源码所在路径。
  4. BUILT_DIR路径,默认值为BUILT_DIR = . ,即当前源码所在路径。

配置好Makefile文件即可编译,使用以下命令进行编译:

make clean
make shud

更多可选的编译参数:

  • make all - 清理旧的编译文件,并编译SHUD和SHUD_omp文件。
  • make help - 展示编译的帮助信息。
  • make shud - 编译shud模型的可执行文件。
  • make shud_omp - 编译OpenMP支持的并行版本,编译后的可执行文件名为shud_omp。

3.1.2.1 OpenMP(可选)

当需要并行计算时,需要安装OpenMP。以下为Mac OS上安装的参数:

brew install llvm clang
brew install libomp
compile flags for OpenMP: 
  -Xpreprocessor -fopenmp -lomp
Library/Include paths:
  -L/usr/local/opt/libomp/lib 
  -I/usr/local/opt/libomp/include

3.1.2.2 运行 SHUD 可执行文件

运行SHUD模型,请在命令行内执行一下命令:

./shud <projectname>

Default screen print of shud and command tips 命令的参数模式:

./shud [-0] [-p projectfile] [-c calibrationfile] [-o output] [-n Num_Threads] <project_name>
  • -0 空计算。读取输入文件并写出结果文件,但是内部没有计算。仅仅用于验证输入文件格式的可靠性和可用性。

  • <project name> 项目名称。所有的项目文件必须以此开头。

  • [-p projectfile] 指定的项目配置文件,文件中包含了所有输入文件的地址。

  • [-c calibrationfile] 指定的参数校准文件(.cfg.calib), 默认的文件为input/xx/xx.cfg.calib。

  • [-o output_folder] 制定的结果文件输出地址。默认的输出地址为output/projname.out

  • [-n Num_Threads] 并行计算的线程数量。此参数仅在并行计算打开时可用。

shud 程序正常运行时,屏幕输出应该如下图所示: Screen print of SHUD running.

3.1.3 rSHUD

rSHUD是一个标准R包,但尚未向CRAN提交,因此需要通过github以源码方式安装。在R环境中执行以下命令即可自动安装rSHUD。 关于rSHUD的内容,.

如果您的电脑上尚未在R环境中安装devtools,请先安装。

install.packages('devtools')

然后通过devtools从github上安装rSHUD。

devtools::install_github('SHUD-System/rSHUD')

rSHUD依赖的R工具包包括:

  • Rcpp
  • reshape2
  • ggplot2
  • gridExtra
  • grid
  • fields
  • xts
  • hydroGOF
  • zoo
  • raster (>= 2.1.0)
  • sp
  • rgeos
  • RTriangle
  • rgdal (>= 1.1.0)
  • proj4
  • abind
  • utils
  • lubridate
  • geometry
  • methods
  • ncdf4
  • GGally
  • doParallel

3.2 模型输入数据

本节从原始数据要求到输入文件内涵,向用户详细介绍模型输入数据。

3.2.1 原始数据

类型 数据名 必要性 要求(格式、步长、数量) 说明
空间数据 流域边界 矢量 连续、封闭、唯一多边形;
空间数据 河流网络 矢量 河流具有方向性;河流无向下分叉;
空间数据 高程 0 单位:米
空间数据 土地利用分类 矢量/栅格 \(n_{lc}\)个分类
空间数据 土壤分类 矢量/栅格 \(n_s\)个分类
空间数据 气象站点 矢量/栅格 \(n_{fc}\)个站点
参数 土地利用 \(n_{lc}\) 土地利用的水力学参数
参数 土壤 \(n_s\) 土壤性质:Silt百分比、Clay百分比、有机质含量、Bulk Density
时间序列数据 气象驱动 小时~日 降雨、气温、湿度、辐射、风速、气压
时间序列数据 叶面积指数(LAI) 小时~月 每类土地利用一个LAI时间序列
时间序列数据 融化系数 小时~月 一个时间序列
时间序列数据 观测数据 小时~月 通常为径流数据、地下水、蒸散发等

3.2.1.1 数据实例

DWR
南苏丹Pongo流域原始数据实例:流域边界(绿色实线),河流(红色实线)、DEM(背景)和气象站点覆盖范围(黑色虚线)。气象数据站点为GLDAS数据,因此为$ 0.25 x 0.25 $度覆盖范围。
DWR DWR
(a) Percentage of silt (b) Percentage of Clay
DWR DWR
(c) Organic mater (Organic Carbon) (d) Bulk density
南苏丹Pongo流域土壤数据:Silt百分比,Clay百分比,有机物含量,Bulk Density。数据来源:ISRIC_SoilGrids。
DWR
南苏丹Pongo流域土地利用信息。土地利用数据来自USGS GLC数据。
DWR
FLDAS数据的覆盖范围,即气象站点位置信息。
DWR
FLDAS的气象驱动数据。

3.2.2 模型输入文件

输入文件列表

File Category 备注 Header # of column
.mesh sp 三角形单元定义
.att sp 单元属性表
.riv sp 河流
.rivseg sp 三角形单元与河道单元拓扑信息
.calib cfg 物理参数校准文件
.para cfg 模型运行设置参数文件
.ic cfg 模型初始条件
.geol para 地下水层的水力学参数
.soil para 浅层土壤的水力学参数
.lc para 土地利用的水力学和能量参数
.forc tsd 气象驱动数据文件的列表
.csv tsd 气象驱动时间序列数据
.lai tsd 叶面积指数
.obs tsd 用于校准观测时间序列数据
.mf tsd 融雪参数时间序列数据
.rl tsd 计算潜在蒸散发的粗糙度时间序列数据
gis/domain Shapefile 三角形单元的矢量文件 x x
gis/river Shapefile 河段的矢量文件 x x
gis/seg Shapefile 被三角形分割的河段矢量文件 x x
The screenshot of input files for SHUD
The screenshot of input files for SHUD

在文件夹 gisfig 中的文件与SHUD模拟无关,然后在数据分析和可视化中有很大作用,因此保留在输入文件当中。

3.2.3 空间数据

3.2.3.1 .sp.mesh 文件

Example of .sp.mesh file (block 1) Example of .sp.mesh file (block 2)

.mesh 文件中有两个表格,第一个表格为三角形的组成和拓扑关系;第二个表格为组成三角形的节点坐标信息。

  • 表格 1 (三角形单元的几何定义)

  • 表头

值1 值2
行数 ( \(N_{cell}\)) 列数 (\(8\))
含义 取值范围 单位 备注
ID 单元序号 \(i\) 1 ~ \(N_{cell}\) -
Node1 三角形单元第1角点\(i\) 1 ~ \(N_{node}\) -
Node2 三角形单元第2角点 \(i\) 1 ~ \(N_{node}\) -
Node3 三角形单元第3角点 \(i\) 1 ~ \(N_{node}\) -
Nabr1 三角形单元第1邻居序号\(i\) 1 ~ \(N_{cell}\) -
Nabr2 三角形单元第2邻居序号\(i\) 1 ~ \(N_{cell}\) -
Nabr3 三角形单元第3邻居序号\(i\) 1 ~ \(N_{cell}\) -
Zmax 三角形单元z中心点地表高程 -9999 ~ +inf \(m\)
  • Block 2 (三角形的三个角点的信息)

  • 表头:

值1 值2
行数( \(N_{node}\)) 列数 (\(5\))
  • Table
含义 取值范围 单位 备注
ID 点序号 \(i\) 1 ~ \(N_{cell}\) -
X x坐标 \(i\) 1 ~ \(N_{node}\) -
Y y坐标 \(i\) 1 ~ \(N_{node}\) -
AqDepth 含水层厚度(地表至不透水层距离)\(i\) 0 ~ +inf \(m\)
Elevation 点的高程(z坐标) \(i\) -9999 ~ +inf \(m\)

3.2.3.2 .sp.att文件

Example of .sp.att file
Example of .sp.att file
  • 表头
值1 值2
行数( \(N_{cell}\)) 列数 (\(7\))
  • Table
含义 取值范围 单位 备注
ID 单元序号\(i\) 1 ~ \(N_{cell}\) -
SOIL 土壤类型序号 1 ~ \(N_{soil}\) -
GEOL 类型序号 1 ~ \(N_{geol}\) -
LC 类型序号 1 ~ \(N_{lc}\) - \(N_{lc}\) = \(N_{lai}\)
FORC 气象站点序号 1 ~ \(N_{forc}\) -
MF 融雪指数序号 1 ~ \(N_{mf}\) -
BC 边界条件序号 1 ~ \(N_{bc}\) -
SS 源汇序号 1 ~ \(N_{bc}\) -

3.2.3.3 .sp.riv文件

Example of .sp.riv file
Example of .sp.riv file
  • 表头
值1 值2
行数( \(N_{riv}\)) 列数 (\(5\))
含义 取值范围 单位 备注
ID 河段序号\(i\) 1 ~ \(N_{river}\) -
DOWN 下游河段序号 1 ~ \(N_{river}\) - Negative vlaue indicates outlet
Type 河流参数序号 1 ~ \(N_{rivertype}\) -
Slope 河床底坡度 -10 ~ 10 \(m/m\) Height/Length
Length 河段长度 \(i\) 0 ~ inf \(m\)

3.2.3.4 .sp.rivseg文件

Example of .sp.rivseg file
Example of .sp.rivseg file
  • 表头
值1 值2
行数( \(N_{segment}\)) 列数 (\(4\))
  • Table
含义 取值范围 单位 备注
ID 片段序号 \(i\) 1 ~ \(N_{segment}\) -
iRiv 所属河段序号 1 ~ \(N_{river}\) -
iEle 相交单元序号 1 ~ \(N_{cell}\) -
Length 片段长度 \(i\) 0 ~ inf \(m\)

3.2.4 水力学参数

3.2.4.1 .para.soil

土壤表层的水力学参数

3.2.4.2 .para.geol

地下水层的水力学参数

3.2.4.3 .para.lc

与土地利用有关的参数

3.2.5 模型配置文件

3.2.5.1 .cfg.para文件

Example of .cfg.para file
Example of .cfg.para file
  • Table
含义 取值范围 单位 推荐值
VERBOSE 输出冗余信息 - - 0
INIT_MODE 初始条件模式 0~3 - 3 (0=Relief condition, 1=Dry condition, 2=Default guess, 3=Warm start)
ASCII_OUTPUT 是否输出文本格式结果 1/0 - 0
Binary_OUTPUT 是否输出二进制格式结果 1/0 - 1
CRYOSPHERE 冰冻圈模块是否开启 1/0 - 0
SPINUPDAY 模型预热天数 0 ~ inf \(day\) 0
SCR_INTV 屏幕输出间隔 0 ~ \(N_{threads}\) \(min\) 1440
ABSTOL CVODE绝对容差 1e-6 ~ 0.1 - 0.0001
RELTOL CVODE相对容差 1e-6 ~ 0.1 - 0.0001
INIT_SOLVER_STEP 初始迭代步长 - \(min\) 1
MAX_SOLVER_STEP 最大迭代步长 1~60 \(min\) 10
ET_STEP 正散发计算步长 1~360 \(min\) 60
START 模型开始时间(天数) 0 ~ inf \(day\) 0
END 模型结束时间(天数) - \(day\) -
dt_ye_snow 储量输出步长,积雪 0 ~ inf \(min\) 1440
dt_ye_surf 储量输出步长,地表水 0 ~ inf \(min\) 1440
dt_ye_unsat 储量输出步长,未饱和层 0 ~ inf \(min\) 1440
dt_ye_gw 储量输出步长,地下水层 0 ~ inf \(min\) 1440
dt_Qe_surf 单元流量输出步长,地表水 0 ~ inf \(min\) 1440
dt_Qe_sub 单元流量输出步长, 地下水 0 ~ inf \(min\) 1440
dt_qe_et0 单元流量输出步长,截流蒸发 0 ~ inf \(min\) 1440
dt_qe_et1 单元流量输出步长,蒸腾 0 ~ inf \(min\) 1440
dt_qe_et2 单元流量输出步长,蒸发 0 ~ inf \(min\) 1440
dt_qe_etp 单元流量输出步长,潜在蒸散发 0 ~ inf \(min\) 1440
dt_qe_prcp 单元流量输出步长,降雨 0 ~ inf \(min\) 1440
dt_qe_infil 单元流量输出步长,下渗 0 ~ inf \(min\) 1440
dt_qe_rech 单元流量输出步长,地下水补给 0 ~ inf \(min\) 1440
dt_yr_stage 河段储量输出步长,河道 0 ~ inf \(min\) 1440
dt_Qr_down 河段流量输出步长,向下游 0 ~ inf \(min\) 1440
dt_Qr_surf 河段流量输出步长,向坡面 0 ~ inf \(min\) 1440
dt_Qr_sub 河段流量输出步长,向地下水 0 ~ inf \(min\) 1440
dt_Qr_up 河段流量输出步长,向上游 0 ~ inf \(min\) 1440
dt_Lake 湖泊相关变量的输出步长 0 ~ inf \(min\) 1440

注:

  • INIT_MODE:
    • 0 = 含水层完全饱和作为初始条件,
    • 1 = 含水层完全干涸,
    • 2 = 未饱和层水量为含水层40%,饱和层水储量为含水层的40%,
    • 3 = 读取.cfg.ic作为初始条件

3.2.5.2 .cfg.calib文件

Example of .cfg.calib file
Example of .cfg.calib file
  • Table
含义 取值范围 单位 备注
GEOL_KSATH 水平水力传导度,地下水层 ? -
GEOL_KSATV 垂直水力传导度,地下水层 ? -
GEOL_KMACSATH 大孔隙水平水力传导度,地下水层 ? -
GEOL_DMAC 大孔隙深度,地下水层 -
GEOL_THETAS 地下水层孔隙度,地下水层 -
GEOL_THETAR 土壤残留含水量,地下水层 -
GEOL_MACVF 大孔隙面积比,地下水层 -
SOIL_KINF 饱和水力传导度,未饱和层 ? -
SOIL_KMACSATV 大孔隙垂向水力传导度,未饱和层 ? -
SOIL_DINF 下渗深度参数,未饱和层 ? -
SOIL_ALPHA \(\alpha\)值,van Genuchten公式 -
SOIL_BETA \(\beta\)值,van Genuchten公式 -
SOIL_MACHF 大孔隙面积比,未饱和层 -
LC_VEGFRAC 植被覆盖度 -
LC_ALBEDO 反照率 -
LC_ROUGH 地表曼宁粗糙度 -
LC_SOILDGD 土壤劣化系数 -
LC_IMPAF 不透水面积比 -
LC_ISMAX 最大截流系数 -
AQ_DEPTH+ 含水层深度 \(m\)
TS_PRCP 降雨 -
TS_SFCTMP+ 温度 \(C\)
ET_ETP 潜在蒸散发 -
ET_IC 冠层截流 -
ET_TR 植被蒸腾 -
ET_SOIL 直接蒸发 -
RIV_ROUGH 河道曼宁粗糙度 -
RIV_KH 河床水力传导度 -
RIV_DPTH+ 河道深度 \(m\)
RIV_WDTH+ 河道宽度 \(m\)
RIV_SINU 河道绵延度 -
RIV_CWR 谢才(Chezy)公式系数\(C_{wr}\) -
RIV_BSLOPE+ 河床坡度 \(m/m\)
IC_GW+ 地下水水位初始条件 \(m\)
IC_RIV+ 河道水量初始条件 \(m\)

注:

  • 表格中+代表该值的含义为该值与模型参数相加得到新参数。即\(\theta ^* = \theta + \zeta\), \(\theta\)为参数, \(\zeta\)为校准调节值。
  • 除+以外的所有值,都是与参数相乘关系。即\(\theta ^* = \theta \times \zeta\)

3.2.5.3 .cfg.ic文件

Example of .cfg.ic file
Example of .cfg.ic file
  • 表1 (三角单元初始条件)

  • 表头

值1 值2
行数( \(N_{cell}\)) 列数 (\(6\))
  • Table
含义 取值范围 单位 备注
ID 单元序号\(i\) 1 ~ \(N_{cell}\) -
Canopy 冠层截流 \(i\) 0 ~ inf \(m\)
Snow 积雪 \(i\) 0 ~ inf \(m\)
Surface 地表水 \(i\) 0 ~ inf \(m\)
Unsat 未饱和层水 \(i\) 0 ~ inf \(m\)
GW 地下水 \(i\) 0 ~ inf \(m\)
  • 表2 (河段初始条件)

  • 表头:

值1 值2
行数( \(N_{riv}\)) 列数 (\(2\))
  • Table
含义 取值范围 单位 备注
ID 河段序号 \(i\) 1 ~ \(N_{riv}\) -
Stage 河段水位 \(i\) 0 ~ inf \(m\)
  • 表3 (湖泊初始条件)

  • 表头:

值1 值2
行数( \(N_{lake}\)) 列数 (\(2\))
  • Table
含义 取值范围 单位 备注
ID 湖泊序号 \(i\) 1 ~ \(N_{lake}\) -
Stage 湖泊水位 \(i\) 0 ~ inf \(m\)

3.2.6 时间序列数据

3.2.6.1 .tsd.forc文件

Example of .tsd.forc file
Example of .tsd.forc file
  • Line 1: 气象站点数量 | 开始日期 (YYYYMMDD)
  • Line 2: 驱动数据文件所在文件夹
  • Line 3~N: 时间序列数据文件名
Example of .csv forcing file
Example of .csv forcing file
  • 表头:
值1 值2
( \(0\)) 列数 (\(6\))
  • Table
含义 取值范围 单位 备注
Day 日数 0 ~ \(N_{day}\) \(day\)
PRCP 降雨 0 ~ 500 \(mm/day\) 降雨强度超过阈值模型会警告,但不妨碍运行
TEMP 气温 -100 ~ 70 \(C\)
RH 相对湿度 0 ~ 1 \(-\)
wind 风速 0 ~ inf \(m/s\) 非负
Rn 太阳净辐射 ? \(W/m^2\) 最大值为太阳常数
Press 气压 ? \(Pa\) 可缺省

3.2.6.2 .tsd.lai文件

Example of .tsd.lai file
Example of .tsd.lai file
  • 表头:
值1 值2 值3
天数 ( \(N_{time}\)) 列数 (\(N_{lc}\)) 开始日期 (YYYYMMDD)
  • Table
含义 取值范围 单位 备注
第一列 时间 0 ~ \(N_{time}\) \(day\)
第2列 叶面积指数1 0 ~ inf \(m^2/m^2\)
第i列 叶面积指数\(i-1\) 0 ~ inf \(m^2/m^2\)

3.2.6.3 .tsd.mf文件

Example of .tsd.mf file
Example of .tsd.mf file

融雪度日模型(Degree-Day Model)的系数。

  • 表头:
值1 值2 值3
天数 ( \(N_{time}\)) 列数 (\(N_{mf}\)) 开始日期 (YYYYMMDD)
  • Table
含义 取值范围 单位 备注
第一列 时间 0 ~ \(N_{time}\) \(day\)
第2列 融雪因子1 0 ~ inf -
第i列 融雪因子\(i-1\) 0 ~ inf -

3.2.6.4 .tsd.obs文件

Example of .tsd.obs file
Example of .tsd.obs file
  • 表头:
值1 值2 值3
天数 ( \(N_{time}\)) 列数 (\(N_{obs}\)) 开始日期 (YYYYMMDD)
  • Table
含义 取值范围 单位 备注
第一列 时间 0 ~ \(N_{time}\) \(day\)
第2列 观测1 ? ?
第i列 观测 \(i-1\) ? ?

3.3 结果输出文件

3.3.1 文件名格式

为了方便结果的分类和用户理解,作者规则化定义了输出文件的文件名;用户通过文件名即可读懂文件所包含的变量类型、所属单元和读写格式。输出文件的文件名格式:

[项目名].[识别码].[格式]

或者:

[Project Name].[Identifier].[format]

  • [项目名][Project Name]) 是用户给模拟定的名称,通常为流域的名字简写。例如:黄河流域,可以写成hh或者huanghe.

  • [格式][format])可以使用csv或者dat. csv为规则文本格式,可使用任意文本编辑器或者WPS表格打开。dat二进制格式,文件占用硬盘空间较小

  • [识别码][Identifier]) 定义了输出结果内涵。识别码基本结构为 [模型单元][变量类型][变量名][Model cell][Variable Type][Variable Name])。 [Model cell]ele (element,模型三角单元), riv (river,河段) or lak (lake,湖泊). 变量类型包括 y, v and q,分别代表状态变量 (单位\(L\)), 单位流量或速度(单位\(L^3 L^{-2} T^{-1}\))和体积流量(in \(L^3 T^{-1}\))——\(L\)为长度单位,\(T\)为时间单位。

下表为SHUD输出文件列表:

识别码 单元 类型 变量名 含义 单位
.eleyic. ele y ic 储量,冠层截流 \(m\)
.eleysnow. ele y snow 储量,雪水当量 \(m\)
.eleysurf. ele y surf 储量,地表水 \(m\)
.eleyunsat. ele y unsat 储量,未饱和层 \(m\)
.eleygw. ele y gw 储量,地下水水位 \(m\)
.elevetp. ele v etp 通量,潜在蒸散发 \(\frac{m^3}{m^2 d}\)
.eleveta. ele v eta 通量,实际蒸散发 \(\frac{m^3}{m^2 d}\)
.elevetic. ele v etic 通量,冠层截流蒸发 \(\frac{m^3}{m^2 d}\)
.elevettr. ele v ettr 通量,植被蒸腾 \(\frac{m^3}{m^2 d}\)
.elevetev. ele v etev 通量,直接蒸发 \(\frac{m^3}{m^2 d}\)
.elevprcp. ele v prcp 通量,降雨 \(\frac{m^3}{m^2 d}\)
.elevnetprcp. ele v netprcp 通量,净降雨量 \(\frac{m^3}{m^2 d}\)
.elevinfil. ele v infil 通量,下渗 \(\frac{m^3}{m^2 d}\)
.elevexfil. ele v infil 通量,出渗(反下渗) \(\frac{m^3}{m^2 d}\)
.elevrech. ele v rech 通量,地下水补给 \(\frac{m^3}{m^2 d}\)
.eleqsurf. ele q surf 流量,地表径流(坡面流) \(m^3/d\)
.eleqsub. ele q sub 流量, 地下水水平流动 \(m^3/d\)
.rivystage. riv y stage 储量,河段水位 \(m\)
.rivqup. riv q up 流量,向上游 \(m^3/d\)
.rivqdown. riv q down 流量,向下游 \(m^3/d\)
.rivqsurf. riv q surf 净流量,坡面向外 \(m^3/d\)
.rivqsurf1. riv q surf 流量,方向1 \(m^3/d\)
.rivqsurf2. riv q surf 流量,方向2 \(m^3/d\)
.rivqsurf3. riv q surf 流量,方向3 \(m^3/d\)
.rivqsub. riv q sub 净流量,地下水向外 \(m^3/d\)
.rivqsub1. riv q sub 流量,方向1 \(m^3/d\)
.rivqsub2. riv q sub 流量,方向2 \(m^3/d\)
.rivqsub3. riv q sub 流量,方向3 \(m^3/d\)

3.3.2 文本格式输出文件(.csv)

N - 输出文件列数(不含时间列) m - 步长总数 StartTime - 开始时间 (YYYYMMDD or YYYYMMDD.hhmmss)

N StartTime
\(T_1\) \(v_{1 \cdot 1}\) \(v_{1 \cdot 2}\) \(v_{1 \cdot N}\)
\(T_2\) \(v_{2 \cdot 1}\) \(v_{2 \cdot 2}\) \(v_{2 \cdot N}\)
\(T_3\) \(v_{3 \cdot 1}\) \(v_{3 \cdot 2}\) \(v_{3 \cdot N}\)
\(T_{m}\) \(v_{m \cdot 1}\) \(v_{m \cdot 2}\) \(v_{m \cdot N}\)

3.3.3 二进制输出文件(.dat)

二进制文件输出的结果核心矩阵与文本格式一致,但是使用了不同的存储结构。二进制文件格式紧凑,读写效率高,占用较小磁盘空间,因此作为SHUD模型的默认输出方式。 二进制文件格式不仅可被各类计算机语言(C/C++, Basic, Fortran,……)快速读取,并且也能在R、python、 Matlab、 Mathmetica等交互式语言中高效读取,如果读取有障碍,请联系模型作者获取帮助。

表格中,重要的变量是数据列数\(N_x\)\(N_x\)代表输出变量的列数,即该变量的输出计算单元数。 \(0 < N_x \le N_{max}\), 例如模型有\(N_{cell}\)个三角形单元,则.eleygw的地下水输出文件中\(N_{max} = N_{cell}\);模型的河流径流结果中\(N_{max} = N_{riv}\)

二进制输出文件现有v1.0, v2.0和v2.1三种格式,三种输出结果的格式见下表:

v 1.0 (\(N_x = N_{max}\)

内容 字节数 类型 示例
数据列数 8 double 100
起始时间 8 double 20220101.120000
数据矩阵 8 double 1440, 1,2,…; 2880,2,4,…

v 2.0(\(N_x = N_{max}\)

内容 字节数 类型 示例
文件描述 1024 char groundwater export, bla bla bla
起始时间 8 double 20220101.120000
数据列数 8 double 100
数据矩阵 8 double 1440, 1,2,…; 2880,2,4,…

v 2.1 以上版本(\(0 < N_x \le N_{max}\)

内容 字节数 类型 示例
版本号 8 double 2.0
文件描述 2048 char groundwater export, bla bla bla
起始时间 8 double 20220101.120000
数据列数 8 double 100
索引ID 8 double 1,3,5,…;
数据矩阵 8 double 1440, 1,2,…; 2880,2,4,…
v2.1版本输出文件的逻辑和内存结构
v2.1版本输出文件的逻辑和内存结构

注:double为计算机的双精度数据格式,默认长度为8个比特位

数据为\(m \times N_x\)的二维矩阵。

  • 矩阵中的值全部为8字节double型。
  • \(m\)为结果输出的最终时间长度,\(m=100\)表明模型共输出100个时间步长的数据,则数据矩阵的行数为100。
  • \(N_x\)为输出结果的数量。 在v2.1之后的版本中,模型可灵活输出指定序号的结果,因此在v2.1结果中,结果矩阵的列数并不一定等于模拟单元的数量;例如模拟三角形单元数为1000(\(N_{max} = 1000\)),模型配置中指定值输出偶数序号的结果,则\(N_x = 500\),索引ID则为\(i_{o} = 2,4,6,...,998,1000\)
  • \(i\)为矩阵列序号,取值范围:\(0 \le i \le N_{x}\),第一列(\(i=0\))为时间,\(i=1~N_x\)则为输出结果值。
  • \(j\)为时间序号,取值范围:$0 j m $。
  • \(N_{h}\)为文件头的长度,换算为8字节double型的数量。

数据矩阵的存储方式见下表:

文件存储位置 矩阵行号(\(j=\)) 数据列序号(\(i=\)) 值含义
\(N_{h} + 1\) 1 0 \(T_1\)
\(N_{h} + 2\) 1 1 \(v_{1 \cdot 1}\)
\(N_{h} + 3\) 1 2 \(v_{1 \cdot 2}\)
1
\(N_{h} +N_x\) 1 \(N_x\) \(v_{1 \cdot N_x}\)
\(N_{h} +(N_x+1) * (2-1) + i\) 2 0 \(T_2\)
\(N_{h} +(N_x+1) * (2-1) + i\) 2 1 \(v_{2 \cdot 1}\)
\(N_{h} +(N_x+1) * (2-1) + i\) 2 2 \(v_{2 \cdot 2}\)
\(N_{h} +(N_x+1) * (2-1) + i\) 2
\(N_{h} +(N_x+1) * (2-1) + i\) 2 \(N_x\) \(v_{2 \cdot N_x}\)
\(N_{h} +(N_x+1) * (3-1) + i\) 3 0 \(T_3\)
\(N_{h} +(N_x+1) * (3-1) + i\) 3 1 \(v_{3 \cdot 1}\)
\(N_{h} +(N_x+1) * (3-1) + i\) 3 2 \(v_{3 \cdot 2}\)
\(N_{h} +(N_x+1) * (3-1) + i\) 3
\(N_{h} +(N_x+1) * (3-1) + i\) 3 \(N_x\) \(v_{3 \cdot N_x}\)
\(N_{h} +(N_x+1) * (j-1) + i\) \(j\)
\(N_{h} +(N_x+1) * (j-1) + i\) \(j\)
\(N_{h} +(N_x+1) * (j-1) + i\) \(j\)
\(N_{h} +(N_x+1) * (j-1) + i\) \(j\)
\(N_{h} +(N_x+1) * (m-1) + i\) m 0 \(T_{m}\)
\(N_{h} +(N_x+1) * (m-1) + i\) m 1 \(v_{m \cdot 1}\)
\(N_{h} +(N_x+1) * (m-1) + i\) m 2 \(v_{m \cdot 2}\)
\(N_{h} +(N_x+1) * (m-1) + i\) m
\(N_{h} +(N_x+1) * (m-1) + i\) m \(N_x\) \(v_{m \cdot N_x}\)