Mechanics Input (mech_3D.inp) (January, 2022) The file mech_3D.inp contains all flags, input and output file names, and geomechanics properties for the geomechanical simulation. File mech_3D.inp is free format file. As for TOUGH2, inputs are preceeded by a keyword heading. Comments and blank lines placed at the end of records are ignored, which allows for detailed user descriptions of input parameters. Elastic properties may be specified for a general horizontally isotropic (vertically anisotropic) elastic medium. However, equations for failure strains assume isotropic material, and any failing material is treated as if it had isotropic elastic properties prior to failure (given by Ep and nu_pp, below). For isotropic material, Ez = Ep, nu_zp = nu_pp, and the Gzp should be Ep / 2 (1+nu_pp) to be consistent with them. order record 1: 'x' 'z' 'y' parameters record 1: nx ny nz Number of mesh blocks in x, y and z directions. record 2: ftol stol ftol - Relative tolerance for solution of mechanics force balance equations (ftol), e.g., 1e-6 . stol - Maximum change in porosity (in one time step) before reducing time step (stol), e.g., 5.e-3. If failure is not expected, stol can be set to a fairly small value (1.e-4) and used to control time stepping. If hydraulic fracturing is expected, a fairly large value (e.g., 5.e-3) should be used. record 3: initstressfile File name for initial stresses. If the file name is given as "smoothedlocalloading" or "uniformstressgradient", no file is read, and the initial stresses are calculated internally. These two cases take additional inputs; record 3a (smoothedlocalloading): initstressratio(1) initstressratio(2) initstressratio(3) horzstressrotation Ratios of sigma_x'x', sigma_y'y', and sigma_zz to computed vertical stress, and rotation of principal stress coordinates x' and y' (rotation of x' from x towards y, in degrees). Vertical stress is computed from local loading by rock mass computed from grain density (flow.inp rocks.1 record, characters 11-20), porosity, and internally computed fluid density (at P, T conditions), and from vertical stress specified at boundary normal to z. Assumes that vertical (zz) stress is set on one surface normal to z. Loading stress at each level is smoothed before integration to next level. When smoothed local loading stresses are computed, they are subsequently written to a file named smoothedlocalloading. record 3a (uniformstressgradient): refstress(1) refstress(2) refstress(3) record 3b (uniformstressgradient): str_grad(1,1) str_grad(1,2) str_grad(1,3) record 3c (uniformstressgradient): str_grad(2,1) str_grad(2,2) str_grad(2,3) record 3d (uniformstressgradient): str_grad(3,1) str_grad(3,2) str_grad(3,3) Diagonal stresses sigma_xx, sigma_yy, sigma_zz (Pascal) at (0,0,0) and gradient of diagonal stresses (Pa/m). Indexing is such that (str_grad(1,1), str_grad(1,2), str_grad(1,3)) = d (sigma_xx, sigma_yy, sigma_zz)/ dx. Compressive stresses are considered negative. record 4: iBC(1,1) iBC(2,1) Mechanics boundary conditions at x(1), x(nx+1). Choices are 0 no displacement normal to boundary 1 no displacement on boundary 2 normal traction on boundary constant to be read in 3 traction on boundary constant value to be read in 4 normal traction on boundary held at initial normal stress in adjacent cell 5 traction on boundary held at initial stress in adjacent cell record 4a (iBC(1,1)= 2 case): surftractn(1,1,1) x traction on x(1) surface. record 4a (iBC(1,1)= 3 case): surftractn(1,1,1) surftractn(2,1,1) surftractn(3,1,1) x,y,z tractions on x(1) surface. Record 4a present when iBC(1,1)= 2 or 3. record 4b (iBC(2,1)= 2 case): surftractn(1,2,1) x traction on x(nx+1) surface. record 4b (iBC(2,1)= 3 case): surftractn(1,2,1) surftractn(2,2,1) surftractn(3,2,1) x,y,z tractions on x(nx+1) surface. Record 4b present when iBC(2,1)= 2 or 3. record 5: iBC(1,2) iBC(2,2) Mechanics boundary conditions at y(1), y(ny+1). Choices 0, 2, 3, 4, or 5, as above. record 5a (iBC(1,2)= 2 case): surftractn(2,1,2) y traction on y(1) surface. record 5a (iBC(1,2)= 3 case): surftractn(1,1,2) surftractn(2,1,2) surftractn(3,1,2) x,y,z tractions on y(1) surface. Record 5a present when iBC(1,2)= 2 or 3. record 5b (iBC(2,2)= 2 case): surftractn(2,2,2) y traction on y(ny+1) surface. record 5b (iBC(2,2)= 3 case): surftractn(1,2,2) surftractn(2,2,2) surftractn(3,2,2) x,y,z tractions on y(ny+1) surface. Record 5b present when iBC(2,2)= 2 or 3. record 6: iBC(1,3) iBC(2,3) Mechanics boundary conditions z(1), z(nz+1). Choices 0, 2, 3, 4, or 5, as above. record 6a (iBC(1,3)= 2 case): surftractn(3,1,3) z traction on z(1) surface. record 6a (iBC(1,3)= 3 case): surftractn(1,1,3) surftractn(2,1,3) surftractn(3,1,3) x,y,z tractions on z(1) surface. Record 6a present when iBC(1,3)= 2 or 3. record 6b (iBC(2,3)= 2 case): surftractn(3,2,3) z traction on z(nz+1) surface. record 6b (iBC(2,3)= 3 case): surftractn(1,2,3) surftractn(2,2,3) surftractn(3,2,3) x,y,z tractions on z(nz+1) surface. Record 6b present when iBC(2,3)= 2 or 3. record 7: 0 0 record 8: 0 0 record 9: 0 0 record 10: BCvol_thresh Threshold for zeroing mechanics porosity changes (e.g., 1.e50): elements with (MESH file) ELEME field volume greater than or equal to BCvol_thresh are not given porosity changes due to mechanics. (Zeroing of porosity changes due to chemistry depends separately on whether ELEME volume is greater than 1020 m3.) petsc_mnp (optional field) record 1: m n p number of processes to assign in x, y and z directions. When program is compiled for MPI parallel processing and run under mpirun, used to override program’s default domain sub-division, otherwise ignored. When not ignored, product m*n*p must match number of MPI processes specified for mpirun. materials record 1: nummat number of materials, i.e., number of rock types. Must be the same number of rock types as in flow.inp. Mechanical properties of each material are given on nummaterlines consecutive lines. record 2: materlinecontrols three 0 or 1's. The first controls whether (1) or not (0) a line of temperature-stress bilinear strain parameters follows the first four lines of material properties for each rock type. The second controls whether or not a line of porosity dependent moduli parameters follows first four (or five) lines of material properties for each rock type. The third controls whether or not a line of volumetric effective stress permeability parameters is read. The number of lines of material properties for each rock type nummaterlines is either 4, 5, 6, or 7, depending on the number of 1's in materlinecontrols (000, 100, 010, 110, 001, 101, 011, or 111). records 3, 3+ nummaterlines , ..., 3 + i*nummaterlines, ...: mname young_par young_perp shear poisson_par poisson_perp angle five character material name (rock type) as in flow.inp ROCKS block, Youngs modulus parallel; Ep (in plane normal to z (anisotropy) direction) (Pascal), Youngs modulus perpendicular; Ez (in anisotropy direction) (Pascal), shear modulus; Gzp for shearing between parallel and z directions, Poisson ratio parallel; nu_pp between x and y directions, Poisson ratio perpendicular; nu_zp ratio of (-)parallel displacement to vertical displacement for an applied vertical stress, (angle of anisotropy from vertical); not used, vertical anisotropic assumed. records 4, 4+ nummaterlines , ..., 4+i x nummaterlines, ...: biot cohesion tensile_stren MC_friction bulk_dens coeff_plast dilation_ang thermal_dilat Biot coefficient, e.g., 1 - K_dr/K_grain . (Value = -1 treated specially, for water elements in which thermal stresses from rock matrix expansion are not present.) Mohr-Coulomb cohesion (Pascal), e.g., 0, tensile strength (Pascal), e.g., 0, Mohr-Coulomb friction angle (degrees), e.g., 30, bulk density (kg/m3) (either this, or flow.inp grain density and porosity is used for mechanics, according to idensitysource below) coefficient of plastic porosity- not used, dilation angle (dilation on M.C. failure) (degrees) thermal dilation coefficient (1/ ?C). records 5, 5 + nummaterlines , ..., 5 + i*nummaterlines, ...: imatrixpermporotype amatrixpermparam bmatrixpermparam cmatrixpermparam ifracpermporotype afracpermparam bfracpermparam index of matrix porosity permeability relationship for rock type (0-6). Types 1-6 as in TOUGHREACT V3.0-OMP. Index 0 defaults to type specified in chemical.inp Zppr record ippzon value, when chemistry is run, and to none otherwise. a_param of matrix porosity permeability relation. b_param of matrix porosity permeability relation. c_param of matrix porosity permeability relation. index of fracture porosity permeability relationship for rock type (0, 7 or 8). 0 neglects permeability changes due to failure porosity. a_param of fracture porosity permeability relation, for ifracpermporotype=7, fraction of initial porosity assumed due to fractures (generally a small number), for ifracpermporotype=8, assumed fracture spacing (m). b_param of fracture porosity permeability relation, for ifracpermporotype=7 or 8, fraction of initial permeability assumed due to fractures (generally close to 1). records 6, 6+ records 6, 6+ nummaterlines , ..., 6 + i*nummaterlines, ...: idensitysource idissolv_poro idensitysource =0 for bulk density from flow.inp grain density and porosity and initial fluid density; =1 for bulk density from here; =2 for bulk density from here and grain density reset from bulk density and flow.inp porosity, assuming a fluid density of 1000 kg/m3 ; =3 for bulk density from here and grain density reset from bulk density and flow.inp porosity, assuming a fluid density of 1000 kg/m3 idissolv_poro = 0, no imposed dissolution strains, = 1, sigma_eff*(T-T_0 )*coef , stress temperature product strains imposed. (idisolv_poro=1 may need quite small time steps for stability.) records 7, 7+ nummaterlines , ..., 6 + i*nummaterlines, ... (when materlinecontrl(1:1)=1): c_bilin T_bilin K_bilin d_bilin bilinear strain parameters (1-4) for this rock type. records 8, 8+ nummaterlines , ..., 7 + i*nummaterlines, … (when materlinecontrl(2:2)=1): phi_crit a_k b_k c_k d_k a_sh b_sh c_sh d_sh porosity dependent moduli parameters for this rock type. Bulk modulus is given by K * ( a_k +b_k (phi -phi_crit )^(1/2) +c_k (phi -phi_crit )+d_k (phi -phi_crit )^2 ) / ( a_k+b_k (phi_0 -phi_crit )^(1/2) +c_k (phi -phi_crit ) + d_k (phi_0 -phi_crit )) and shear moduls by mu * ( a_sh +b_sh (phi -phi_crit )^(1/2) +c_sh (phi -phi_crit ) +d_sh (phi-phi_crit )^2 ) / (a_sh+b_sh (phi_0- phi_crit )^(1/2) +c_sh (phi -phi_crit ) + d_sh (phi_0- phi_crit )^2 ), where K and mu are the starting bulk and shear moduli. records 9, 9+ nummaterlines , ..., 8 +i*nummaterlines, ... (when materlinecontrl(3:3)=1): log10(sigma_k0) dlog10(sigma_efftrans)/dlog10(T) dlog10(k)/dlog10(sigma_vol)|low dlog10(k)/dlog10(sigma_v)|high volumetric effective stress dependent permeability parameters, after Watanabe, et al. (2017). Transistion volumetric effective stress is eigma_efftrans = log10(sigma_k0) + log10(T) * dlog10(sigma_efftrans)/dlog10(T). For sigma_effvol <= sigma_efftrans have delta log10(k) = dlog10(k)/dlog10(sigma_vol)|low & delta log10(sigma_effvol). For sigma_effvol > sigma_efftrans have delta log10(k) = dlog10(k)/dlog10(sigma_vol)|high delta log10(sigma_effvol). Watanabe et al.'s values for these parameters are 12.8, -4.3, -0.9, and -3.8. Here, T is temperature Celsius, k is permeability in m^2, and, with respect to these parameters, volumetric stressess are in MPa, with compressive stress positive, both contrary to the rest of this file. When using this option, for materials for which it is not to be applied, set values to 0 0 0 0. Generally, when using this option, one wants imatrixpermporotype=0 for the materials it is applied to. bilinearstrain (optional field) record 1: c_bilin T_bilin K_bilin d_bilin tau_smoo For use with rock types with idisolv_poro=1 (above). Sets parameters for temperature -volumetric effective stress product dependent strain of form: epsilon_d = sigma_eff [ ( T - T_bilin) c_bilin - 1/(3 K_bilin) ] - epsilon_d0 where epsilon_d0 = sigma_eff(0) [ ( T(0) - T_bilin) c_bilin - 1/(3 K_bilin) ] when indicator trace(sigma_eff ) == (T-T_bilin) is increasing, that is, when volumetric effective stress magnitude trace (sigma_eff) or temperature T increases. When this indicator decreases, strain epsilon_d reverses, but with only d_bilin as much change as indicated above (0 <= d_bilin <= 1). These strains are smoothed over time scale ?_smoo for stability. When materials field nmaterlines < 5, all five bilinearstrain entries here are used. When materials nmaterlines =5, only fifth bilinearstrain parameter (tau_smoo) is used from here, other values are taken from fifth line of material properties for each rock type. bilinearstrain parameters have no effect on rock types with idissolv_poro = 0. gridmaterials (optional field) record 1: ngridmater number (<=14) of mechanics properties to take from geomechanprop.inp file instead of from mech_3D.inp materials block. record 2: igridmater_1 igridmater_2 ... igridmater_ ngridmater indices (between 1 and 14) of material properties to take from geomechanprop.inp instead of from mech_3D.inp materials block. (igridmater=1: young_par, igridmater=14: thermal_dilat.) (Properties 15 to 22 are always set by rock type, and properties 23-26 either by rock type or by the bilinearstrain field parameters.) mopm record 1: mopmech(1) mopmech(2) mopmech(3) mopmech(4) mopmech5) numMINCMe iverbose ForPoroC chm_rel_p chm_rel_T maxiter_me ibiot_eff_stress ifailureporo poromix pelastmix mopmech(1) = 0 no mechanics; = 1 mechanics, =2 read mechanics input files and stop. mopmech(2) -- not used mopmech(3) -- not used mopmech(4) = 0 write total stress, biot*pressure, temp, total strain only at elements requested by elgm (below). (RMECH.out files produced with mopmech(4)=0 are insufficient for use as an initial stresses file.) = 1 write total stress, biot*pressure, temp, total strain, average displacement at every element, to file RMECH.out every time step = 2 write total stress, biot*pressure, temp, total strain, average displacement at every element, to file RMECH.out when save files are updated = 3 write total stress, biot*pressure, temp, total strain at every element to file RMECH.out each time step, overwriting old file = 4 write total stress, biot*pressure, temp, total strain at every element to file RMECH.out after first time step and stop = 5 write initstressratio(1:3) * smoothed local loading stresses to file smoothedlocalloading and stop Cases 0, 3, 4, & 5 are the most practical, as RMECH.out files from them either grow slowly (case 0), or remain a constant size. When ibiot_eff_stress = 0 (below), pressure is written to file RMECH.out in place of biot*pressure. mopmech(5) not used numMINCMe number of MINC levels (0 no MINC) iverbose = 0 quiet > 0 diagnostic information written to standard output mopmech(8) =0, default, =1, on restarts, zero displacements, strains, failure strains, and dissolution strains, (leaves stresses unchanged) =2, on restarts, same as 1, plus, restores any cohesion lost on failure chm_rel_p -- not used chm_rel_T -- not used maxiter_me -- not used ibiot_eff_stress = 0 use unweighted pressure in effective stress for Mohr-Coulomb and tensile failure criteria (following Dempsey et al., 2014, and Jaeger, et al., 2007), = 1 weight pressure by Biot coefficient in computing effective stress for Mohr-Coulomb and tensile failure criteria (Zoback, 2014, undated lecture notes). ifailureporo = 0 neglect changes in porosity due to material failure when computing pressure, = 1 compute pressure taking into account changes in porosity due to material failure. poromix fraction of change in failure porosity (e.g., 0.25) to incorporate immediately in pressure calculations. poromix < 1 increases stability under hydraulic fracturing conditions. Failure porosity changes get averaged in in 1/poromix time steps (in a `leaky bucket' auto-regressive average). When hydraulic fracturing is not expected to occur, poromix = 1 is fine. pelastmix fraction of change in poroelastic porosity (e.g., 0.25) to incorporate immediately in pressure calculations. pelastmix < 1 increases stability under hydraulic fracturing conditions. Poroelastic porosity changes get averaged in in 1/pelastmix time steps (in a `leaky bucket' auto-regressive average). When hydraulic fracturing is not expected to occur, pelastmix = 1 is fine. elgm (optional field) record 1: numprelem number of mechanics elements to write Monitor_CNNNN files for, where CNNNN corresponds to flow.inp element name. record 2: monelem_mechname (1) ELEME name (in MESH file) of element to write mechanics monitor information, ... record 1+numprelem: monelem_mechname (numprelem) ELEME name (in MESH file) of element to write mechanics monitor information, permlimit (optional field) record 1: ixpermlimit iypermlimit izpermlimit indices of element to center perm limit function on. record 2: radiuspermlimit permlim_atradius radius (m) of uniform maximum permeability limit zone. Permeability limit increases as r for r > radiuspermlimit. maximum allowed permeability (m2) in uniform maximum permeability zone. record 3; nexemptrocktype number of rock types exempt from permeability limit. record 4; rtype name of first rock type to be exempt (not present for nexemptrocktype=0). ... record 3+nexemptrocktype; rtype name of last rock type to be exempt from permeability limit. vertaverage (optional field) record 1; iz1vertavg iz2vertavg z indices of elements to average stresses and strains between for output file stress_strain_vertavg.tec done (optional field) end of mech_3D.inp file. Mechanics initial stresses file. The geomechanics initial stresses file name is read from mech_3D.inp file parameters field record 3. The initial stresses file has free format. The default format of the file is; Title time 1 sigma_xx (1,1) sigma_xx (2,1) sigma_xx (3,1) sigma_xx (4,1) sigma_xx (5,1) sigma_xx (6,1) sigma_xx (7,1) sigma_xx (8,1) sigma_yy (1,1) sigma_yy (2,1) sigma_yy (3,1) sigma_yy (4,1) sigma_yy (5,1) sigma_yy (6,1) sigma_yy (7,1) sigma_yy (8,1) sigma_zz (1,1) sigma_zz (2,1) sigma_zz (3,1) sigma_zz (4,1) sigma_zz (5,1) sigma_zz (6,1) sigma_zz (7,1) sigma_zz (8,1) sigma_yz (1,1) sigma_yz (2,1) sigma_yz (3,1) sigma_yz (4,1) sigma_yz (5,1) sigma_yz (6,1) sigma_yz (7,1) sigma_yz (8,1) sigma_xz (1,1) sigma_xz (2,1) sigma_xz (3,1) sigma_xz (4,1) sigma_xz (5,1) sigma_xz (6,1) sigma_xz (7,1) sigma_xz (8,1) sigma_xy (1,1) sigma_xy (2,1) sigma_xy (3,1) sigma_xy (4,1) sigma_xy (5,1) sigma_xy (6,1) sigma_xy (7,1) sigma_xy (8,1) 2 sigma_xx (1,2) sigma_xx (2,2) sigma_xx (3,2) sigma_xx (4,2) sigma_xx (5,2) sigma_xx (6,2) sigma_xx (7,2) sigma_xx (8,2) sigma_yy (1,2) sigma_yy (2,2) sigma_yy (3,2) sigma_yy (4,2) sigma_yy (5,2) sigma_yy (6,2) sigma_yy (7,2) sigma_yy (8,2) sigma_zz (1,2) sigma_zz (2,2) sigma_zz (3,2) sigma_zz (4,2) sigma_zz (5,2) sigma_zz (6,2) sigma_zz (7,2) sigma_zz (8,2) sigma_yz (1,2) sigma_yz (2,2) sigma_yz (3,2) sigma_yz (4,2) sigma_yz (5,2) sigma_yz (6,2) sigma_yz (7,2) sigma_yz (8,2) sigma_xz (1,2) sigma_xz (2,2) sigma_xz (3,2) sigma_xz (4,2) sigma_xz (5,2) sigma_xz (6,2) sigma_xz (7,2) sigma_xz (8,2) sigma_xy (1,2) sigma_xy (2,2) sigma_xy (3,2) sigma_xy (4,2) sigma_xy (5,2) sigma_xy (6,2) sigma_xy (7,2) sigma_xy (8,2) ... nelem sigma_xx (1,nelem) sigma_xx (2,nelem) sigma_xx (3,nelem) sigma_xx (4,nelem) sigma_xx (5,nelem) sigma_xx (6,nelem) sigma_xx (7,nelem) sigma_xx (8,nelem) sigma_yy (1,nelem) sigma_yy (2,nelem) sigma_yy (3,nelem) sigma_yy (4,nelem) sigma_yy (5,nelem) sigma_yy (6,nelem) sigma_yy (7,nelem) sigma_yy (8,nelem) sigma_zz (1,nelem) sigma_zz (2,nelem) sigma_zz (3,nelem) sigma_zz (4,nelem) sigma_zz (5,nelem) sigma_zz (6,nelem) sigma_zz (7,nelem) sigma_zz (8,nelem) sigma_yz (1,nelem) sigma_yz (2,nelem) sigma_yz (3,nelem) sigma_yz (4,nelem) sigma_yz (5,nelem) sigma_yz (6,nelem) sigma_yz (7,nelem) sigma_yz (8,nelem) sigma_xz (1,nelem) sigma_xz (2,nelem) sigma_xz (3,nelem) sigma_xz (4,nelem) sigma_xz (5,nelem) sigma_xz (6,nelem) sigma_xz (7,nelem) sigma_xz (8,nelem) sigma_xy (1,nelem) sigma_xy (2,nelem) sigma_xy (3,nelem) sigma_xy (4,nelem) sigma_xy (5,nelem) sigma_xy (6,nelem) sigma_xy (7,nelem) sigma_xy (8,nelem) where each stress component is specified at the eight Gauss points in each element. There are nelem lines of (48) stresses, and nelem = nx * ny * nz. Time is in seconds and stresses are in Pascals. Extra items at the ends of lines are ignored. Files with title Geomechanics Restart are treated specially, assuming all information written to a geomechan.restart file is present. For an element with corners at (x(i),y(j),z(k)), (x(i+1),y(j),z(k)), (x(i),y(j),z(k+1)), (x(i+1),y(j),z(k+1)), (x(i),y(j+1),z(k)), (x(i+1),y(j+1),z(k)), (x(i),y(j+1),z(k+1)), (x(i+1),y(j+1),z(k+1)), Gauss points are 1/[2 sqrt(3)] grid spacings out, in each dimension, from the cell center, towards corners (x(i),y(j),z(k)), (x(i+1),y(j),z(k)), (x(i+1),y(j+1),z(k)), (x(i),y(j+1),z(k)), (x(i),y(j),z(k+1)), (x(i+1),y(j),z(k+1)), (x(i+1),y(j+1),z(k+1)), (x(i),y(j+1),z(k+1)), respectively. Stresses at elements are in MESH file element order. Initial stresses file names smoothedlocalloading and uniformstressgradient in mech_3D.inp are treated specially and do not require files of those names. Instead, they invoke reading additional lines of input from the mech_3D.inp file and either calculate stresses internally from local loading stresses, or impose a uniform gradient diagonal stress field, respectively. More details are given in mech_3D.inp file format above.