-
Notifications
You must be signed in to change notification settings - Fork 148
MODFLOW NetCDF Format
This guide describes the supported formats of MODFLOW 6 NetCDF files. As NetCDF support in MODFLOW 6 is an ongoing development effort, aspects of these formats related specifically to MODFLOW 6 (naming conventions, internal attributes, etc.) are subject to change.
Extended MODFLOW 6 can read and write version 4 (NetCDF-4) NetCDF files. Two different formats are supported for each (input and output) file type: a structured format and a UGRID based unstructured LAYERED MESH format. Each format aims to be cf-conventions compliant.
All MODFLOW 6 NetCDF files correspond to a model. Output files typically contain arrays for the model dependent variable. Input files can contain package arrays that are supported as inputs from NetCDF files. These variables are indicated with red highlighting in MODFLOW 6 - Description of Input and Output (mf6io.pdf). Currently, supported inputs include GRIDDATA block arrays and PERIOD block arrays for READASARRAYS and READARRAYGRID based packages.
MODFLOW 6 structured NetCDF files support GWF, GWT and GWE DIS based simulations. These files minimally define the x, y, z, and time dimensions for data variables. The time dimension corresponds to the total number of simulation steps. x, y, and z correspond to ncol, nrow, and nlay of the model discretization. GRIDDATA arrays are dimensioned (z, y, x), while PERIOD arrays are dimensioned (time, z, y, x).
MODFLOW 6 UGRID Layered Mesh NetCDF files support GWF, GWT and GWE DIS and DISV based simulations. Each grid based array variable in these file types is split into a set of layer variables. NPF K input data, for example, would be defined in 10 layer data variables when the grid has 10 layers. These files also differ from the structured variants because they contain a set of variables whose purpose is to geometrically define the grid. This is adventageous for displaying the grid itself as well as data associated with the grid in applications that support reading the mesh, such as QGIS.
UGRID layered mesh files minimally define the dimensions time, nmesh_node, nmesh_face, max_nmesh_face_nodes, x and y. The time dimension corresponds to the total number of simulation steps. nmesh_node and nmesh_face correspond to number of grid vertices and ncpl, respectively. max_nmesh_face_nodes coorespondes to the maximum number of vertices in a cell, 4 when modflow_grid is STRUCTURED. x and y are defined for input data that is associated with the row or column dimensions of a layer, for example DELC and DELR. GRIDDATA arrays are dimensioned (nmesh_face), while PERIOD arrays are dimensioned (time, nmesh_face).
MODFLOW 6 NetCDF input files must define MODFLOW 6 specific attributes to be processed correctly. These attributes are defined at both the file (dataset) and variable scope and vary according to the file format. Any variable that is not annotated with the expected attributes will be ignored and as such MODFLOW 6 NetCDF input files can contain data not intended for MODFLOW 6.
There are no specific requirements related to variable naming other than from the conventions themselves.
Two file scoped attributes are required for structured NetCDF inputs: modflow_grid and modflow_model. modflow_grid describes the MODFLOW 6 discretization type and modflow_model describes the model type and name. For structured inputs, modflow_grid should always be set to STRUCTURED. The modflow_model attribute should be set to MODFLOW model type: model name, for example: GWF6: GWF-VSC01-BND
Data variables intended to be read as MODFLOW 6 package input must be annotated correctly with MODFLOW 6 input attributes. Variables without the expected annotation will be ignored. The modflow_input attribute is always required while the modflow_iaux attribute is required for auxiliary variables. Use of these attributes are described in more detail below.
File scoped attributes modflow_grid and modflow_model are required as described for the structured input. DIS based simulation input should set modflow_grid to STRUCTURED, while DISV based simulation input should set modflow_grid to VERTEX. An additional mesh attribute is required for UGRID based NetCDF input files. UGRID Layered Mesh files should set this attribute to LAYERED.
Data variables intended to be read as MODFLOW 6 package input must be annotated correctly with MODFLOW 6 input attributes. Variables without the expected annotation will be ignored. The modflow_input attribute is always required while the modflow_iaux attribute is required for auxiliary variables. Use of these attributes are described in more detail below. In addition, UGRID layered mesh input parameters must define the layer attribute and set it to the interger value of the layer that the data is assoicated with.
The modflow_input variable scoped attribute is always required and corresponds to the variable model, package and input tag name. The model name should match the name in the modflow_model file scoped attribute. This name should also match a model name in the mfsim.nam models block. The package name is the MODFLOW 6 package name or type, described next. Base packages should use the package type as name, e.g. NPF or DISV. Stress packages should use the package name in the model name file packages block, or the FloPy generated name when the packages block does not define a name.
The modflow_iaux variable is an additional attribute required only for auxiliary array variables. MODFLOW 6 auxiliary variables are configured with the AUXILIARY keyword in supported options blocks. This keyword is associated with a list of auxiliary names, for example TEMPERATURE followed by SALINITY. modflow_iaux should be set to the position in the auxiliary list for each input aux array. Here, any input data associated with the TEMPERATURE auxiliary should set modflow_iaux to 1 while any SALINITY array should set modflow_iaux to 2.
Here are ascii dumps using the ncdump netcdf command line utility. The dumps represent netcdf binary file dataset headers. Dataset dimensions are listed in the top section followed by a section of variables. Each variable definition shows it's dimensions and attributes, including MODFLOW 6 specific attributes. Dataset scoped attributes are in the last section of the dump.
As the dump from the following model netcdf input file has twenty layers, the number of variables in a URID based represention of the same input data would be more than 20 times higher, which would include grid variables. The size of the files would be more comparable, however, as the amount of user input data would be the same.
netcdf gwf_henrynr01 {
dimensions:
bnd = 2 ;
time = 1250 ;
z = 20 ;
y = 1 ;
x = 40 ;
variables:
double time(time) ;
time:calendar = "standard" ;
time:units = "days since 1970-01-01T00:00:00" ;
time:axis = "T" ;
time:standard_name = "time" ;
time:long_name = "time" ;
double z(z) ;
z:units = "layer" ;
z:long_name = "layer number" ;
double y(y) ;
y:units = "m" ;
y:axis = "Y" ;
y:standard_name = "projection_y_coordinate" ;
y:long_name = "Northing" ;
y:bounds = "y_bnds" ;
double y_bnds(y, bnd) ;
double x(x) ;
x:units = "m" ;
x:axis = "X" ;
x:standard_name = "projection_x_coordinate" ;
x:long_name = "Easting" ;
x:bounds = "x_bnds" ;
double x_bnds(x, bnd) ;
double wel-1_q(time, z, y, x) ;
wel-1_q:_FillValue = 3.e+30 ;
wel-1_q:units = "m" ;
wel-1_q:long_name = "well rate" ;
wel-1_q:modflow_input = "GWF_HENRYNR01/WEL-1/Q" ;
double wel-1_concentration(time, z, y, x) ;
wel-1_concentration:_FillValue = 3.e+30 ;
wel-1_concentration:units = "m" ;
wel-1_concentration:long_name = "well auxiliary variable iaux concentration" ;
wel-1_concentration:modflow_input = "GWF_HENRYNR01/WEL-1/AUX" ;
wel-1_concentration:modflow_iaux = 1 ;
double drn-1_elev(time, z, y, x) ;
drn-1_elev:_FillValue = 3.e+30 ;
drn-1_elev:units = "m" ;
drn-1_elev:long_name = "drain elevation" ;
drn-1_elev:modflow_input = "GWF_HENRYNR01/DRN-1/ELEV" ;
double drn-1_cond(time, z, y, x) ;
drn-1_cond:_FillValue = 3.e+30 ;
drn-1_cond:units = "m" ;
drn-1_cond:long_name = "drain conductance" ;
drn-1_cond:modflow_input = "GWF_HENRYNR01/DRN-1/COND" ;
double drn-1_concentration(time, z, y, x) ;
drn-1_concentration:_FillValue = 3.e+30 ;
drn-1_concentration:units = "m" ;
drn-1_concentration:long_name = "drain auxiliary variable iaux concentration" ;
drn-1_concentration:modflow_input = "GWF_HENRYNR01/DRN-1/AUX" ;
drn-1_concentration:modflow_iaux = 1 ;
double ghb-1_bhead(time, z, y, x) ;
ghb-1_bhead:_FillValue = 3.e+30 ;
ghb-1_bhead:units = "m" ;
ghb-1_bhead:long_name = "boundary head" ;
ghb-1_bhead:modflow_input = "GWF_HENRYNR01/GHB-1/BHEAD" ;
double ghb-1_cond(time, z, y, x) ;
ghb-1_cond:_FillValue = 3.e+30 ;
ghb-1_cond:units = "m" ;
ghb-1_cond:long_name = "boundary conductance" ;
ghb-1_cond:modflow_input = "GWF_HENRYNR01/GHB-1/COND" ;
double ghb-1_concentration(time, z, y, x) ;
ghb-1_concentration:_FillValue = 3.e+30 ;
ghb-1_concentration:units = "m" ;
ghb-1_concentration:long_name = "general-head boundary auxiliary variable iaux concentration" ;
ghb-1_concentration:modflow_input = "GWF_HENRYNR01/GHB-1/AUX" ;
ghb-1_concentration:modflow_iaux = 1 ;
double ghb-1_density(time, z, y, x) ;
ghb-1_density:_FillValue = 3.e+30 ;
ghb-1_density:units = "m" ;
ghb-1_density:long_name = "general-head boundary auxiliary variable iaux density" ;
ghb-1_density:modflow_input = "GWF_HENRYNR01/GHB-1/AUX" ;
ghb-1_density:modflow_iaux = 2 ;
// global attributes:
:title = "GWF_HENRYNR01 hydraulic head array input" ;
:source = "MODFLOW 6 6.7.0.dev2 (preliminary) 05/12/2025" ;
:modflow_grid = "STRUCTURED" ;
:modflow_model = "GWF6: GWF_HENRYNR01" ;
:history = "first created 2025/8/2 21:1:38.17" ;
:Conventions = "CF-1.11" ;
}
Here is an example UGRID format input file for a GWT model with 1 layer:
netcdf gwt_dsp01a {
dimensions:
time = 200 ;
nmesh_node = 202 ;
nmesh_face = 100 ;
max_nmesh_face_nodes = 4 ;
x = 100 ;
y = 1 ;
variables:
double time(time) ;
time:calendar = "standard" ;
time:units = "days since 2041-01-01T00:00:00-05:00" ;
time:axis = "T" ;
time:standard_name = "time" ;
time:long_name = "time" ;
int mesh ;
mesh:cf_role = "mesh_topology" ;
mesh:long_name = "2D mesh topology" ;
mesh:topology_dimension = 2 ;
mesh:face_dimension = "nmesh_face" ;
mesh:node_coordinates = "mesh_node_x mesh_node_y" ;
mesh:face_coordinates = "mesh_face_x mesh_face_y" ;
mesh:face_node_connectivity = "mesh_face_nodes" ;
double mesh_node_x(nmesh_node) ;
mesh_node_x:units = "m" ;
mesh_node_x:standard_name = "projection_x_coordinate" ;
mesh_node_x:long_name = "Easting" ;
double mesh_node_y(nmesh_node) ;
mesh_node_y:units = "m" ;
mesh_node_y:standard_name = "projection_y_coordinate" ;
mesh_node_y:long_name = "Northing" ;
double mesh_face_x(nmesh_face) ;
mesh_face_x:units = "m" ;
mesh_face_x:standard_name = "projection_x_coordinate" ;
mesh_face_x:long_name = "Easting" ;
mesh_face_x:bounds = "mesh_face_xbnds" ;
double mesh_face_xbnds(nmesh_face, max_nmesh_face_nodes) ;
double mesh_face_y(nmesh_face) ;
mesh_face_y:units = "m" ;
mesh_face_y:standard_name = "projection_y_coordinate" ;
mesh_face_y:long_name = "Northing" ;
mesh_face_y:bounds = "mesh_face_ybnds" ;
double mesh_face_ybnds(nmesh_face, max_nmesh_face_nodes) ;
int mesh_face_nodes(nmesh_face, max_nmesh_face_nodes) ;
mesh_face_nodes:cf_role = "face_node_connectivity" ;
mesh_face_nodes:long_name = "Vertices bounding cell (counterclockwise)" ;
mesh_face_nodes:_FillValue = -2147483647 ;
mesh_face_nodes:start_index = 1 ;
double dis_delr(x) ;
dis_delr:_FillValue = 9.96920996838687e+36 ;
dis_delr:long_name = "spacing along a row" ;
dis_delr:modflow_input = "GWT_DSP01A/DIS/DELR" ;
double dis_delc(y) ;
dis_delc:_FillValue = 9.96920996838687e+36 ;
dis_delc:long_name = "spacing along a column" ;
dis_delc:modflow_input = "GWT_DSP01A/DIS/DELC" ;
double dis_top(nmesh_face) ;
dis_top:_FillValue = 9.96920996838687e+36 ;
dis_top:long_name = "cell top elevation" ;
dis_top:modflow_input = "GWT_DSP01A/DIS/TOP" ;
double dis_botm_l1(nmesh_face) ;
dis_botm_l1:_FillValue = 9.96920996838687e+36 ;
dis_botm_l1:long_name = "cell bottom elevation layer=1" ;
dis_botm_l1:modflow_input = "GWT_DSP01A/DIS/BOTM" ;
dis_botm_l1:layer = 1 ;
int dis_idomain_l1(nmesh_face) ;
dis_idomain_l1:_FillValue = -2147483647 ;
dis_idomain_l1:long_name = "idomain existence array layer=1" ;
dis_idomain_l1:modflow_input = "GWT_DSP01A/DIS/IDOMAIN" ;
dis_idomain_l1:layer = 1 ;
double ic_strt_l1(nmesh_face) ;
ic_strt_l1:_FillValue = 9.96920996838687e+36 ;
ic_strt_l1:long_name = "starting concentration layer=1" ;
ic_strt_l1:modflow_input = "GWT_DSP01A/IC/STRT" ;
ic_strt_l1:layer = 1 ;
double dsp_diffc_l1(nmesh_face) ;
dsp_diffc_l1:_FillValue = 9.96920996838687e+36 ;
dsp_diffc_l1:long_name = "effective molecular diffusion coefficient layer=1" ;
dsp_diffc_l1:modflow_input = "GWT_DSP01A/DSP/DIFFC" ;
dsp_diffc_l1:layer = 1 ;
double dsp_alh_l1(nmesh_face) ;
dsp_alh_l1:_FillValue = 9.96920996838687e+36 ;
dsp_alh_l1:long_name = "longitudinal dispersivity in horizontal direction layer=1" ;
dsp_alh_l1:modflow_input = "GWT_DSP01A/DSP/ALH" ;
dsp_alh_l1:layer = 1 ;
double dsp_alv_l1(nmesh_face) ;
dsp_alv_l1:_FillValue = 9.96920996838687e+36 ;
dsp_alv_l1:long_name = "longitudinal dispersivity in vertical direction layer=1" ;
dsp_alv_l1:modflow_input = "GWT_DSP01A/DSP/ALV" ;
dsp_alv_l1:layer = 1 ;
double dsp_ath1_l1(nmesh_face) ;
dsp_ath1_l1:_FillValue = 9.96920996838687e+36 ;
dsp_ath1_l1:long_name = "transverse dispersivity in horizontal direction layer=1" ;
dsp_ath1_l1:modflow_input = "GWT_DSP01A/DSP/ATH1" ;
dsp_ath1_l1:layer = 1 ;
double dsp_atv_l1(nmesh_face) ;
dsp_atv_l1:_FillValue = 9.96920996838687e+36 ;
dsp_atv_l1:long_name = "transverse dispersivity when flow is in vertical direction layer=1" ;
dsp_atv_l1:modflow_input = "GWT_DSP01A/DSP/ATV" ;
dsp_atv_l1:layer = 1 ;
// global attributes:
:title = "GWT_DSP01A concentration array input" ;
:source = "MODFLOW 6 6.7.0.dev2 (preliminary) 05/12/2025" ;
:modflow_grid = "STRUCTURED" ;
:mesh = "LAYERED" ;
:modflow_model = "GWT6: GWT_DSP01A" ;
:history = "first created 2025/8/2 21:16:41.722" ;
:Conventions = "CF-1.11 UGRID-1.0" ;
}