Step 18: Choose boundary conditions
A boundary condition must be specified for:
Regional models: every node around the margin of the grid, and these nodes must be ordered in a counterclockwise direction.
Global models: every node on the footwall of
a "subduction zone" thrust fault. Boundary conditions are required
for these nodes to represent possible effects of subducting slabs of
lithosphere, which are not included in the domain of any Shells model.
Subduction zones are defined as faults of less than a critical dip angle
(variable SUBDIP of Shells.f90,
currently 19° ).
Fortunately, a previous step (#16) explained a method for obtaining an
ordered list of boundary nodes, in which each line already includes:
sequence-number, node-number, latitude, longitude.
Next, you must edit this list to create a boundary conditions (.bcs) file. Use any word processor that can save plain ASCII text with no formatting.
The first line of .bcs file should be a descriptive title, which will be copied into output (.out) files, and then copied from there into maps of the solution. A good choice of title is a brief reference to a publication which presents a plate-tectonic model, if one exists.
There are as many following lines as there are boundary nodes. (The node at the beginning of the circuit is not listed again at the end.) In each line, insert a boundary condition index code between the node-number and the (latitude, longitude) coordinate pair:
|
Boundary-condition index code (ICOND) |
Velocity boundary condition |
Traction boundary condition |
Additional data following the code (and before the latitude & longitude) |
|
-1 |
(free) |
no shear tractions; |
(none) |
|
0 |
(free) |
no shear tractions; |
(none) |
|
1 |
velocity component along specified azimuth is fixed; orthogonal component is free |
may be estimated from consistent nodal reaction force anomalies after solution |
velocity magnitude (in SI, this would be in m/s), then azimuth (in degrees clockwise from local North) |
|
2 |
velocity is completely fixed, to be along the specified azimuth at the specified rate |
may be estimated from consistent nodal reaction force anomalies after solution |
velocity magnitude (in SI, this would be in m/s), then azimuth (in degrees clockwise from local North) |
|
3 [N.B. Only tested with subduction zones in global models. Also, now deprecated.] |
velocity component in the direction of subduction is fixed; orthogonal component is free. Note that subduction direction is the direction of subducting plate velocity in the current reference frame (specified by 2-byte input parameter PLTREF [see previous page]), according to the PB2002 rigid-plate model. Therefore, if you use any type-3 boundary conditions, your stress and strain rate results will NOT be invariant to changes in PLTREF! Another issue is that type-3 boundary conditions cannot be used for nodes which are part of plate PLTREF, because its subduction direction is undefined. |
may be estimated from consistent nodal reaction force anomalies after solution |
(none) |
|
4 |
automatic computation of velocity of subducting plate, according to rigid-plate model PB2002, in the current velocity reference frame. (Note: Use code 4 only for footwall nodes of subduction zones!) |
may be estimated from consistent nodal reaction force anomalies after solution |
(none) |
|
5 |
automatic computation of velocity of selected plate, according to rigid-plate model PB2002, in the current velocity reference frame |
may be estimated from consistent nodal reaction force anomalies after solution |
2-byte plate code, enclosed in quotation marks (e.g., “AF”, or “NA”, or “PA”) as given on map in Step 17. |
Notice that users who want to apply fixed boundary velocities according to model PB2002 of Bird [2003] do not have to do any Euler computations; Shells will do them for you; just apply boundary-condition codes 4 and/or 5. (Code 4 would be used for footwalls of subduction zones, and code 5 would be used for side boundaries.)
However, if you want to use different fixed velocities around model edges, you might want to compute these from Euler poles different than those in PB2002. You would then use boundary-condition code(s) 1 and/or 2, and you would have to supply the velocity magnitudes and azimuths for each boundary node. In order to compute the velocities and azimuths needed for type-1 and type-2 boundary conditions, you will need some tools. Individual velocities for a few nodes can be computed with my utility program OmegaXR. If there are many boundary nodes with these index codes, it may be tedious to type in all their coordinates individually. In that case, you can use my program BCs_Tool to process them in batches. This program requires you to temporarily enter boundary condition index code 3 for all the nodes you want to compute that are associated with a single Euler pole, then specify that Euler pole as the program runs. The .bcs file will be re-written with all code-3 indexes converted to code-2, with the necessary velocities and azimuths supplied. If you wish, you can then manually change code-2 to code-1.
If your model is regional and has a side boundary, you need to be aware of some subtleties of the topology. First, if you have placed fault elements along the perimeter of your model, then the boundary nodes are the ones outside the fault, and they belong the adjacent plate, not the one whose volume you are modeling. Therefore, assign them the velocities of the neighboring plate(s). Second, wherever there are faults along the boundary, there is the possibility that they connect in triple-junctions to additional faults which are outside the perimeter of the model. That is, at every end of every boundary fault element, there is a chance of a change in the name and Euler pole of the neighboring plate. For this reason, the grid will have two distinct nodes on the outside of the boundary at each point where two boundary faults meet. Both require boundary conditions. The first one listed belongs the boundary fault which comes first as you go counterclockwise around the boundary.