Step 32: Choose initial values of weighting parameters L₀, A₀, m, ξ

The heart of the NeoKinema algorithm is an (iterated) weighted-least-squares solution for the long-term-average horizontal velocities of all nodes.
The global-misfit function (or, “objective function”) which is minimized is a sum of four types of misfit:

· Continuum (i.e., long-term-average {non-elastic; permanent} strain-rates different from zero in unfaulted elements)

· Stress (i.e., azimuths of most-compressive principal strain-rate different from interpolated stress directions; see last Step)

· Offset-rate (i.e., differences between geologic target offset-rates of faults, and model rates)

· Geodetic (i.e., long-term-average benchmark velocities different from GPS rates {after correction for model seismicity rates})

Equation (1) of Bird [2009] gives the equation for this objective function.
Its shows how each misfit rate is first made dimensionless, by dividing it by an appropriate datum uncertainty (standard deviation, sigma, σ).
It also shows how all dimensionless misfits are squared (so their signs become unimportant, and only their sizes matter).
Finally, it shows how these 4 kinds of misfit terms differ in dimensionality:

· Continuum and Stress misfits are area-integrals over the total surface area of modeled lithosphere.

· Offset-rate misfits are line-integrals along the traces of the relevant active faults.

· Geodetic misfits are point values at benchmarks.

Therefore, the objective function MUST include two dimensional weighting factors:

· L₀, a reference length in m

· A₀, a reference area in m²

whose inverses multiply the {Offset-rate misfit term} and the {Continuum + Stress misfit terms}, respectively,
in order to make all terms in the objective function dimensionless and comparable to each other.

Loosely speaking, L₀ is the “length of fault trace whose misfit has the same weight as one geodetic benchmark.”

Loosely speaking, A₀ is the “area of unfaulted lithosphere whose misfits have the same weights as one geodetic benchmark.”

Another important parameter is m (Greek mu, or mu_ in Fortran), which is the uncertainty (standard error)
in the prior assumption of zero long-term-average {non-elastic; permanent} strain-rate for all unfaulted lithosphere.
It has units of strain-rate; that is, inverse seconds (/s or s^-1), and represents an “average” or “typical” (e.g., RMS)
strain-rate that cannot be eliminated.  Some physical causes contributing to m  include terminations of faults
without connections to other faults; triple-junction instability at 3-fault intersections, and also true distributed
deformation caused by strain-rate-strengthening viscous response of the lithosphere to stress (even at low T).
Other contributions to m probably come from model imperfections, such as an incomplete map of active faults.

[Experts on NeoKinema will understand that m is allowed to have different (pre-assigned) values in different elements,
and that F-E-grid editor OrbWin can display and adjust these values.  However, beginners in NeoKinema modeling
are advised to use only a single, uniform value of m, which they specify in the NeoKinema parameter file p*.nki.]

Finally, one more obscure weighting factor is ξ (Greek xi, or xi_ in Fortran), which is also in units of strain-rate (/s).
It is a small “quantum” of strain-rate which we are willing to accept as a numerical error, as we try to match all the data
(and prior constraints) that NeoKinema is struggling with.  Theoretically, we might like to have ξ go to zero;
however, experience shows that numerical artifacts will appear (due to the F-E approximation, and also due
to numerical errors in the solutions of linear systems) if ξ is forced to be too small.  Typically, we want ξ < m .

All 4 of these weighting parameters are required to be positive.

The initial comments in the NeoKinema_v5.3.f90.txt file (which are worth reading, by the way…)
include some suggestions for values of these 4 parameters that you might try in your p*.nki file:

    2.00E4        L0 = length of fault trace whose offset rate gets unit weight (in m)
    8.00E8        A0 = area of continuum whose stiffness & isotropy get unit weight (in m**2)
    40            refinements (for nonlinearity) allowed in this solution (try 20-40)
    5.0E-16       mu_ = scalar measure of typical anelastic strain rates in continuum (/s)
    1.0E-17       xi_ = small strain-rate increment, in /s  (Not TOO small!  Suggestions: 1.0E-17 to 1.0E-20 /s.)

These values are based on results of Bird [2009] in modeling neotectonics of the western United States.
Perhaps they may also be good starting values for other places on Earth?
However, if you model the lithospheres of other planets (or moons) where strain-rates are thought to be much higher or lower
than on Earth, it might be appropriate to scale m and ξ in proportion.