Transient vadose zone flow and Tc-99 transport#
This tutorial illustrates Amanzi simulation of transient vadose zone flow and non-reactive (tracer) transport in the context of Tc-99 migration at the DOE Hanford BC Cribs and Trenches site. The tutorial involves a two-dimensional heterogeneous three-layer subsurface system with time- and space-varying infiltration at the ground surface.
BC Cribs and Trenches site#
The two-dimensional geometry and salient features of the Hanford BC Cribs and Trenches site, as represented in this simplified tutorial, are defined in the following schematic diagram:
General recharge at the site (\(U\)) increased following development of the Hanford facility in 1956, and Cribs B-17 and B-18 represent localized sources of water infiltration and Tc-99 contamination (\(C_{Tc-99}\)) to the subsurface:
ID |
\(U\) (mm/yr) |
\(C_{Tc-99}\) (mol/m3) |
|---|---|---|
General site, pre-1956 |
3.5 |
N/A |
General site, post-1956 |
47 |
N/A |
B-17, Jan 1956 |
8025 |
1.88e-6 |
B-18, Feb-Mar 1956 |
10439 |
2.27e-6 |
The three facies underlying the site have the following material properties:
Material property |
Facies 1 |
Facies 2 |
Facies 3 |
|---|---|---|---|
Porosity |
0.4082 |
0.2206 |
0.2340 |
Particle density (kg/m3) |
2720.0 |
2720.0 |
2720.0 |
Horizontal permeability (m2) |
1.9976e-12 |
6.9365e-11 |
2.0706e-09 |
Vertical permeability (m2) |
1.9976e-13 |
6.9365e-12 |
2.0706e-09 |
Capillary pressure model |
van Genuchten |
van Genuchten |
Brooks-Corey |
Relative permeability model |
Mualem |
Mualem |
Mualem |
\(\alpha\) (Pa-1) |
1.9467e-04 |
2.0260e-03 |
2.0674e-03 |
\(S_r\) |
0.0 |
0.0 |
0.0 |
van Genuchten \(m\) |
0.2294 |
0.2136 |
N/A |
Brooks-Corey \(\lambda\) |
N/A |
N/A |
0.3006 |
Amanzi input specifications#
A working Amanzi input file for this tutorial simulation is listed in a
separate section below.
Amanzi input file structure and XML element names are intended to be reasonably
self-explanatory; however, selective commentary is provided here for further guidance.
Currently Amanzi requires the SI units indicated in the model_description/units
element. The corresponding units of density, viscosity, pressure,
hydraulic conductivity, and permeability are thus
kg/m3, Pa \(\cdot\) s, Pa, m/s, and m2, respectively.
Simulation outputs also use these dimensional units.
The use of named constants in the definitions element, e.g.,
<constant name="pre_1956_recharge" type="area_mass_flux" value="1.1071E-7"/>
allows the user to use a fixed name throughout the input file, but change its
numerical value in only one location. Here the time- and space-varying infiltration
parameters are input as defined constants. time_macro name="Macro 1" defines the
output times for visualization in seconds; for example, the last time is
6.33834396E10 seconds or year 2008.5 for a 365.25 day year.
Macro 2 is defined but not used in this
example. Under process_kernels the chemistry kernel is turned off because the
simulation involves only non-reactive transport. In the phases
element, equation-of-state (eos) computations are turned off in favor of direct
specification of fluid density and viscosity. Under the mesh element
a three-dimensional grid is specified, however, only one cell is specified
for the \(y\)-coordinate (ny = "1") making the simulation effectively
2D in the \((x,z)\) coordinates. The grid spacing is 0.5 meters
horizontally and 0.42 meters in the vertical dimension. In the
initial_conditions element, a hydrostatic profile is specified through a
linear pressure profile with slope/gradient = \(- \rho g\) = -9789 Pa/m.
The reference point for \(P = P_{atm}\) = 101325 Pa is set to
\(z\) = 0.5 meters rather than the water table elevation of
\(z\) = 0.0 meters to be compatible with a prior benchmarking comparison to
a STOMP code model. The boundary conditions
at the top surface are defined in the problem statement as a volumetric fluxes
[m3/m2d = m/d]. Amanzi currently requires boundary
conditions of this type to be specified as mass fluxes,
\(\rho U\) [kg /m2s]. The inward_mass_flux values
in the XML input have these units.
Amanzi execution#
Amanzi is executed from the command-line using this or analogous command for the user’s specific installation:
mpirun -n 4 /ascem/amanzi/install/current/bin/amanzi --xml_schema=/ascem/amanzi/install/current/bin/amanzi.xsd --xml_file=dvz_3_layer_2d-isv2.xml
Here execution using 4 processor cores is specified by -n 4. Successful completion
is marked by SIMULATION_SUCCESSFUL followed by a timing summary
Amanzi::SIMULATION_SUCCESSFUL
**********************************************************
*** Timing Summary ***
**********************************************************
Amanzi results#
The figure shows saturation dynamics. Saturation grows underneath the left and right cribs during their operational cycles. Note it takes time for two moving plumes to penetrate into the middle soil which leads to saturation increase on the interface between soils.
The next figure shows concentration of Tc-99. The concentration plumes are moving almost vertically with slight latter interaction due to a non-zero pressure gradient in this direction. The speed of the plumes is reduced significantly in the middle soil. The concentration satisfies the maximum principle.
Amanzi XML input file#
<amanzi_input version="1.6-dev" type="unstructured">
<model_description name="BC Cribs">
<comments>Three homogeneous layers in 2D with two cribs discharging Tc-99</comments>
<model_id>dvz_3_layer_2d</model_id>
<author>Vicky Freedman, David Mouton, Erin Barker</author>
<units>
<length_unit>m</length_unit>
<time_unit>s</time_unit>
<mass_unit>kg</mass_unit>
<conc_unit>molar</conc_unit>
</units>
</model_description>
<definitions>
<constants>
<constant name="zero" type="none" value="0.000"/>
<constant name ="start" type="time" value="1956.0;y"/>
<constant name ="B-17_release_begin" type="time" value="1956.0;y"/>
<constant name ="B-17_release_end" type="time" value ="1956.0822;y"/>
<constant name ="B-18_release_begin" type="time" value ="1956.1644;y"/>
<constant name ="B-18_release_end" type="time" value ="1956.3288;y"/>
<numerical_constant name="pre_1956_recharge" type="area_mass_flux" value="1.1071E-7"/>
<numerical_constant name="post_1956_recharge" type="area_mass_flux" value="1.48666E-6"/>
<numerical_constant name="future_recharge" type="area_mass_flux" value="1.48666E-6"/>
</constants>
<macros>
<time_macro name="Macro 1">
<time>6.17266656E10</time>
<time>6.172982136E10</time>
<time>6.173297712E10</time>
<time>6.173613288E10</time>
<time>6.173928864E10</time>
<time>6.17424444E10</time>
<time>6.174560016E10</time>
<time>6.174875592E10</time>
<time>6.175191168E10</time>
<time>6.175506744E10</time>
<time>6.17582232E10</time>
<time>6.176453472E10</time>
<time>6.177084624E10</time>
<time>6.177715776E10</time>
<time>6.178346928E10</time>
<time>6.17897808E10</time>
<time>6.18213384E10</time>
<time>6.1852896E10</time>
<time>6.2010684E10</time>
<time>6.2168472E10</time>
<time>6.232626E10</time>
<time>6.2484048E10</time>
<time>6.2641836E10</time>
<time>6.2799624E10</time>
<time>6.2957412E10</time>
<time>6.31152E10</time>
<time>6.3272988E10</time>
<time>6.3372710016E10</time>
<time>6.33834396E10</time>
</time_macro>
<time_macro name = "Macro 2">
<time>6.17266656E10</time>
<time>6.33834396E10</time>
</time_macro>
<cycle_macro name = "Every_1000_timesteps">
<start>0</start>
<timestep_interval>1000</timestep_interval>
<stop>-1 </stop>
</cycle_macro>
</macros>
</definitions>
<!-- if flow,transport,chemistry are missing, this implies that state == off -->
<process_kernels>
<comments>This is a proposed comment field for process_kernels</comments>
<flow state = "on" model = "richards"/>
<transport state = "on" algorithm = "explicit first-order" sub_cycling = "on"/>
<chemistry state = "off" engine = "none"/>
</process_kernels>
<phases>
<comments>Eliminated term "Uniform" from viscosity and density input. Designed for additional phases to be included.</comments>
<liquid_phase name = "water">
<eos model="constant" viscosity="1.002e-3" density="998.2"/>
<molar_mass>0.018</molar_mass>
<dissolved_components>
<solutes>
<solute coefficient_of_diffusion="1.0e-9">Tc-99</solute>
</solutes>
</dissolved_components>
</liquid_phase>
</phases>
<execution_controls>
<verbosity level="high" />
<execution_control_defaults init_dt="10" max_dt="20 y" reduction_factor="0.8"
increase_factor="1.25" mode="transient" method="bdf1"/>
<execution_control start="0.0;y" end="1956.0,y" init_dt= "1000" max_dt="500.0 y" reduction_factor="0.8" mode="steady"/>
<execution_control start="B-17_release_begin" />
<execution_control start="B-17_release_end" />
<execution_control start="B-18_release_begin" />
<execution_control start="B-18_release_end" end="3000.0;y" />
</execution_controls>
<!-- There are few different ways to iterate to steady-state -
most use pseudo-time. But what may complicate this is that two different PKs
may have different needs here or in transient.-->
<numerical_controls>
<comments>Numerical controls comments here</comments>
<steady-state_controls>
<comments>Note that this section contained data on timesteps, which was moved into the execution control section.</comments>
<min_iterations>10</min_iterations>
<max_iterations>15</max_iterations>
<max_preconditioner_lag_iterations>30</max_preconditioner_lag_iterations>
<nonlinear_tolerance>1.0e-5</nonlinear_tolerance>
</steady-state_controls>
<transient_controls>
<comments>Proposed comments section.</comments>
<bdf1_integration_method min_iterations="10" max_iterations="15" max_preconditioner_lag_iterations="5" />
</transient_controls>
<linear_solver>
<comments>Proposed comment section.</comments>
<method>gmres</method>
<max_iterations>50</max_iterations>
<tolerance>1.0e-18</tolerance>
<preconditioner name = "hypre_amg">
<hypre_cycle_applications>10</hypre_cycle_applications>
<hypre_smoother_sweeps>3</hypre_smoother_sweeps>
<hypre_tolerance>0.1</hypre_tolerance>
<hypre_strong_threshold>0.4</hypre_strong_threshold>
</preconditioner>
</linear_solver>
</numerical_controls>
<regions>
<box name="All" low_coordinates= "0.0,0.0,0.0" high_coordinates="216.0,1.0,107.52"/>
<box name="Between_Planes_0_and_1" low_coordinates= "0.0,0.0,0.0" high_coordinates="216.0,1.0,39.9"/>
<box name="Between_Planes_1_and_2" low_coordinates= "0.0,0.0,39.9" high_coordinates="216.0,1.0,80.22"/>
<box name="Between_Planes_2_and_3" low_coordinates= "0.0,0.0,80.22" high_coordinates="216.0,1.0,107.52"/>
<box name="Crib_216-B-17" low_coordinates= "74.5,0.0,107.52" high_coordinates="78.5,1.0,107.52"/>
<box name="Crib_216-B-18" low_coordinates= "143.5,0.0,107.52" high_coordinates="147.5,1.0,107.52"/>
<box name="Recharge_Boundary_WestofCribs" low_coordinates= "0.0,0.0,107.52" high_coordinates="74.5,1.0,107.52"/>
<box name="Recharge_Boundary_btwnCribs" low_coordinates= "78.5,0.0,107.52" high_coordinates="143.5,1.0,107.52"/>
<box name="Recharge_Boundary_EastofCribs" low_coordinates= "147.5,0.0,107.52" high_coordinates="216.0,1.0,107.52"/>
<box name="Water Table" low_coordinates= "0.0,0.0,0.0" high_coordinates="216.0,1.0,0.0"/>
</regions>
<mesh framework="mstk"> <!-- default is MSTK for unstructured -->
<comments>Pseudo 2D</comments>
<dimension>3</dimension>
<generate>
<number_of_cells nx="432" ny="1" nz="256"/>
<box low_coordinates="0.0,0.0,0.0" high_coordinates="216.0,1.0,107.52"/>
</generate>
</mesh>
<materials>
<material name="Facies_1">
<comments>Material corresponds to region facies1</comments>
<mechanical_properties>
<porosity model="constant" value="0.4082"/>
<particle_density value="2720.0"/>
</mechanical_properties>
<permeability x="1.9976E-12" y="1.9976E-12" z="1.9976E-13"/>
<cap_pressure model="van_genuchten">
<parameters m="0.2294" alpha="1.9467E-04" sr="0.0"/>
</cap_pressure>
<rel_perm model="mualem">
</rel_perm>
<assigned_regions>Between_Planes_1_and_2</assigned_regions>
</material>
<material name ="Facies_2">
<comments>Material corresponds to region facies2</comments>
<mechanical_properties>
<porosity model="constant" value="0.2206" />
<particle_density value="2720.0"/>
</mechanical_properties>
<permeability x="6.9365E-11" y="6.9365E-11" z="6.9365E-12"/>
<cap_pressure model="van_genuchten">
<parameters m="0.2136" alpha="2.0260E-03" sr="0.0"/>
</cap_pressure>
<rel_perm model="mualem">
</rel_perm>
<assigned_regions>Between_Planes_0_and_1</assigned_regions>
</material>
<material name ="Facies_3">
<comments>Material corresponds to region facies3</comments>
<mechanical_properties>
<porosity model="constant" value="0.234" />
<particle_density value="2720.0"/>
</mechanical_properties>
<permeability x="2.0706E-09" y="2.0706E-09" z="2.0706E-09"/>
<cap_pressure model="brooks_corey">
<parameters lambda="0.3006" alpha="2.0674E-03" sr="0.0"/>
</cap_pressure>
<rel_perm model="mualem">
<optional_krel_smoothing_interval>0.0</optional_krel_smoothing_interval>
</rel_perm>
<assigned_regions>Between_Planes_2_and_3</assigned_regions>
</material>
</materials>
<initial_conditions>
<initial_condition name="Pressure and concentration for entire domain">
<comments>Initial Conditions Comments</comments>
<assigned_regions>All</assigned_regions>
<liquid_phase name = "water">
<liquid_component name = "water">
<linear_pressure value = "101325" reference_coord ="0.0,0.0,0.5" gradient="0.0,0.0,-9793.5192"/>
</liquid_component>
<solute_component name = "Tc-99" value = "0" function="uniform"/>
</liquid_phase>
<temperature formula="300.0"/>
</initial_condition>
</initial_conditions>
<boundary_conditions>
<comments>Format was changed because it is more readable for long time series</comments>
<boundary_condition name = "Recharge at top of the domain">
<comments>Would be nice to have comments section for each boundary_condition but its not clear how it maps to UI</comments>
<assigned_regions>Recharge_Boundary_WestofCribs,Recharge_Boundary_btwnCribs,Recharge_Boundary_EastofCribs</assigned_regions>
<liquid_phase name = "water">
<liquid_component name = "water">
<inward_mass_flux start="0.0" function= "constant" value="pre_1956_recharge"/>
<inward_mass_flux start="1956.0,y" function= "constant" value="post_1956_recharge"/>
<inward_mass_flux start="2012.0,y" function= "constant" value="future_recharge"/>
<inward_mass_flux start="3000.0,y" function= "constant" value="future_recharge"/>
</liquid_component>
<solute_component name = "solute">
<aqueous_conc name = "Tc-99" start="0.0" function= "constant" value="zero"/>
<aqueous_conc name = "Tc-99" start="1956.0,y" function= "constant" value="zero"/>
<aqueous_conc name = "Tc-99" start="2012.0,y" function= "constant" value="zero"/>
<aqueous_conc name = "Tc-99" start="3000.0,y" function= "constant" value="zero"/>
</solute_component>
</liquid_phase>
</boundary_condition>
<boundary_condition name = "B-17">
<assigned_regions>Crib_216-B-17</assigned_regions>
<liquid_phase name = "water">
<liquid_component name = "water">
<inward_mass_flux start="0.0" function="constant" value="pre_1956_recharge"/>
<inward_mass_flux start="B-17_release_begin" function="constant" value="0.00254022"/>
<inward_mass_flux start="B-17_release_end" function="constant" value="post_1956_recharge"/>
<inward_mass_flux start="2012.0,y" function="constant" value="future_recharge"/>
<inward_mass_flux start="3000.0,y" function="constant" value="future_recharge"/>
</liquid_component>
<solute_component name = "solute">
<aqueous_conc name="Tc-99" start="0.0" function="constant" value="0.0"/>
<aqueous_conc name="Tc-99" start="B-17_release_begin" function="constant" value="1.881389E-06"/>
<aqueous_conc name="Tc-99" start="B-17_release_end" function="constant" value="0.0"/>
<aqueous_conc name="Tc-99" start="2012.0,y" function="constant" value="0.0"/>
<aqueous_conc name="Tc-99" start="3000.0,y" function="constant" value="0.0"/>
</solute_component>
</liquid_phase>
</boundary_condition>
<boundary_condition name = "B-18">
<assigned_regions>Crib_216-B-18</assigned_regions>
<liquid_phase name = "water">
<liquid_component name = "water">
<inward_mass_flux start="0.0" function= "constant" value="pre_1956_recharge"/>
<inward_mass_flux start="1956.0,y" function= "constant" value="post_1956_recharge"/>
<inward_mass_flux start="B-18_release_begin" function= "constant" value="0.00330423"/>
<inward_mass_flux start="B-18_release_end" function= "constant" value="post_1956_recharge"/>
<inward_mass_flux start="2012.0,y" function= "constant" value="future_recharge"/>
<inward_mass_flux start="3000.0,y" function= "constant" value="future_recharge"/>
</liquid_component>
<solute_component name = "solute">
<aqueous_conc name = "Tc-99" start="0.0" function= "constant" value="0.0"/>
<aqueous_conc name = "Tc-99" start="1956.0,y" function= "constant" value="0.0"/>
<aqueous_conc name = "Tc-99" start="B-18_release_begin" function= "constant" value="2.266885E-06"/>
<aqueous_conc name = "Tc-99" start="B-18_release_end" function= "constant" value="0.0"/>
<aqueous_conc name = "Tc-99" start="2012.0,y" function= "constant" value="0.0"/>
<aqueous_conc name = "Tc-99" start="3000.0,y" function= "constant" value="0.0"/>
</solute_component>
</liquid_phase>
</boundary_condition>
<boundary_condition name = "Bottom of Domain">
<assigned_regions>Water Table</assigned_regions>
<liquid_phase name = "water">
<liquid_component name = "water">
<uniform_pressure name = "1" start="0.0" function= "constant" value="101325.0"/>
<uniform_pressure name = "2" start="3000.0,y" function= "constant" value="101325.0"/>
</liquid_component>
<solute_component name = "solute">
<aqueous_conc name = "Tc-99" start="0.0" function= "constant" value="0.0"/>
<aqueous_conc name = "Tc-99" start="3000.0,y" function= "constant" value="0.0"/>
</solute_component>
</liquid_phase>
</boundary_condition>
</boundary_conditions>
<output>
<vis>
<base_filename>plot</base_filename>
<num_digits>5</num_digits>
<time_macro>Macro 1</time_macro>
</vis>
<checkpoint>
<base_filename>chk</base_filename>
<num_digits>5</num_digits>
<cycle_macro>Every_1000_timesteps</cycle_macro>
</checkpoint>
</output>
</amanzi_input>