We present a general purpose solver for quadratic programs based on the alternating direction method of multipliers, employing a novel operator splitting technique that requires the solution of a quasi-definite linear system with the same coefficient matrix in each iteration. Our algorithm is very robust, placing no requirements on the problem data such as positive definiteness of the objective function or linear independence of the constraint functions. It is division-free once an initial matrix factorization is carried out, making it suitable for real-time applications in embedded systems. In addition, our technique is the first operator splitting method for quadratic programs able to reliably detect primal and dual infeasible problems from the algorithm iterates. The method also supports factorization caching and warm starting, making it particularly efficient when solving parametrized problems arising in finance, control, and machine learning. Our open-source C implementation OSQP has a small footprint, is library-free, and has been extensively tested on many problem instances from a wide variety of application areas. It is typically ten times faster than competing interior point methods, and sometimes much more when factorization caching or warm start is used.

OSQP: An Operator Splitting Solver for Quadratic Programs

Crust at many divergent plate boundaries forms primarily by the injection of vertical sheet-like dykes, some tens of kilometres long1. Previous models of rifting events indicate either lateral dyke growth away from a feeding source, with propagation rates decreasing as the dyke lengthens2, 3, 4, or magma flowing vertically into dykes from an underlying source5, 6, with the role of topography on the evolution of lateral dykes not clear. Here we show how a recent segmented dyke intrusion in the Barðarbunga volcanic system grew laterally for more than 45 kilometres at a variable rate, with topography influencing the direction of propagation. Barriers at the ends of each segment were overcome by the build-up of pressure in the dyke end; then a new segment formed and dyke lengthening temporarily peaked. The dyke evolution, which occurred primarily over 14 days, was revealed by propagating seismicity, ground deformation mapped by Global Positioning System (GPS), interferometric analysis of satellite radar images (InSAR), and graben formation. The strike of the dyke segments varies from an initially radial direction away from the Barðarbunga caldera, towards alignment with that expected from regional stress at the distal end. A model minimizing the combined strain and gravitational potential energy explains the propagation path. Dyke opening and seismicity focused at the most distal segment at any given time, and were simultaneous with magma source deflation and slow collapse at the Barðarbunga caldera, accompanied by a series of magnitude M > 5 earthquakes. Dyke growth was slowed down by an effusive fissure eruption near the end of the dyke. Lateral dyke growth with segment barrier breaking by pressure build-up in the dyke distal end explains how focused upwelling of magma under central volcanoes is effectively redistributed over long distances to create new upper crust at divergent plate boundaries.

/pdf/segmented-lateral-dyke-growth-in-a-rifting-event-at-5bg334vaou.pdf

Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland

Faults in Earth’s crust accommodate slow relative motion between tectonic plates through either similarly slow slip or fast, seismic-wave-producing rupture events perceived as earthquakes. These types of behaviour are often assumed to be separated in space and to occur on two different types of fault segment: one with stable, rate-strengthening friction and the other with rate-weakening friction that leads to stick-slip. The 2011 Tohoku-Oki earthquake with moment magnitude M_w = 9.0 challenged such assumptions by accumulating its largest seismic slip in the area that had been assumed to be creeping. Here we propose a model in which stable, rate-strengthening behaviour at low slip rates is combined with coseismic weakening due to rapid shear heating of pore fluids, allowing unstable slip to occur in segments that can creep between events. The model parameters are based on laboratory measurements on samples from the fault of the M_w 7.6 1999 Chi-Chi earthquake. The long-term slip behaviour of the model, which we examine using a unique numerical approach that includes all wave effects, reproduces and explains a number of both long-term and coseismic observations—some of them seemingly contradictory—about the faults at which the Tohoku-Oki and Chi-Chi earthquakes occurred, including there being more high-frequency radiation from areas of lower slip, the largest seismic slip in the Tohoku-Oki earthquake having occurred in a potentially creeping segment, the overall pattern of previous events in the area and the complexity of the Tohoku-Oki rupture. The implication that earthquake rupture may break through large portions of creeping segments, which are at present considered to be barriers, requires a re-evaluation of seismic hazard in many areas.

Stable creeping fault segments can become destructive as a result of dynamic weakening

Nonvolcanic tremor is observed in close association with geodetically observed slow-slip events in subduction zones. Accumulating evidence points to these events as members of a family of slow earthquakes that occur as shear slip on the downdip extensions of fault zones in a regime that is transitional between a frictionally locked region above and a freely slipping region below. By virtue of their locations and their properties, slow earthquakes are certain to provide new insights into the behavior of earthquakes and faulting and into the hazard they embody.

Slow Earthquakes and Nonvolcanic Tremor

If model parameterizations of unresolved physics, such as the variety of upper ocean mixing processes, are to hold over the large range of time and space scales of importance to climate, they must be strongly physically based. Observations, theories, and models of oceanic vertical mixing are surveyed. Two distinct regimes are identified: ocean mixing in the boundary layer near the surface under a variety of surface forcing conditions (stabilizing, destabilizing, and wind driven), and mixing in the ocean interior due to internal waves, shear instability, and double diffusion (arising from the different molecular diffusion rates of heat and salt). Mixing schemes commonly applied to the upper ocean are shown not to contain some potentially important boundary layer physics. Therefore a new parameterization of oceanic boundary layer mixing is developed to accommodate some of this physics. It includes a scheme for determining the boundary layer depth h, where the turbulent contribution to the vertical shear of a bulk Richardson number is parameterized. Expressions for diffusivity and nonlocal transport throughout the boundary layer are given. The diffusivity is formulated to agree with similarity theory of turbulence in the surface layer and is subject to the conditions that both it and its vertical gradient match the interior values at h. This nonlocal “K profile parameterization” (KPP) is then verified and compared to alternatives, including its atmospheric counterparts. Its most important feature is shown to be the capability of the boundary layer to penetrate well into a stable thermocline in both convective and wind-driven situations. The diffusivities of the aforementioned three interior mixing processes are modeled as constants, functions of a gradient Richardson number (a measure of the relative importance of stratification to destabilizing shear), and functions of the double-diffusion density ratio, Rρ. Oceanic simulations of convective penetration, wind deepening, and diurnal cycling are used to determine appropriate values for various model parameters as weak functions of vertical resolution. Annual cycle simulations at ocean weather station Papa for 1961 and 1969–1974 are used to test the complete suite of parameterizations. Model and observed temperatures at all depths are shown to agree very well into September, after which systematic advective cooling in the ocean produces expected differences. It is argued that this cooling and a steady salt advection into the model are needed to balance the net annual surface heating and freshwater input. With these advections, good multiyear simulations of temperature and salinity can be achieved. These results and KPP simulations of the diurnal cycle at the Long-Term Upper Ocean Study (LOTUS) site are compared with the results of other models. It is demonstrated that the KPP model exchanges properties between the mixed layer and thermocline in a manner consistent with observations, and at least as well or better than alternatives.

Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization

[1] The mechanics of slow slip events (SSE) in subduction zones remain unresolved. We suggest that SSE nucleate in areas of unstable friction under drained conditions, but as slip accelerates dilatancy reduces pore pressure p quenching instability. Competition between dilatant strengthening and thermal pressurization may control whether slip is slow or fast. We model SSE with 2-D elasticity, rate-state friction, and a dilatancy law where porosity evolves toward steady state ss over distance dc and ss = 0 + ln(v/v0); v is slip speed. We consider two diffusion models. Membrane diffusion (MD) is approximated by −(p − p∞)/tf where p and p∞ are shear zone and remote pore pressure and tf is a characteristic diffusion time. Homogeneous diffusion (HD) accurately models fault-normal flow with diffusivity chyd. For MD, linearized analysis defines a boundary ℰ = 1 − a/b between slow and fast slip, where ℰ ≡ f0e/βb(σ − p∞), f0, a, and b are friction parameters and β is compressibility. When ℰ 1 − a/b slip speeds remain quasi-static. For HD, Ep ≡ eh/(β (σ − p∞)) defines dilatancy efficiency, where h is shear zone thickness and v∞ is plate velocity. SSE are favored by large eh and low effective stress. The ratio Ep to thermal pressurization efficiency scales with 1/(σ − p∞), so high p∞ favors SSE, consistent with seismic observations. For Ep ∼ 10−3 transient slip rates, repeat times, average slip, and stress drops are comparable to field observations. Model updip propagation speeds are comparable to those observed along-strike. Many simulations exhibit slow phases driven by steady downdip slip and faster phases that relax the accumulated stress. Model SSE accommodate only a fraction of plate motion; the remaining deficit must be accommodated during coseismic or postseismic slip.

/pdf/dilatant-strengthening-as-a-mechanism-for-slow-slip-events-tab5w63ltn.pdf

Dilatant strengthening as a mechanism for slow slip events

[1] Episodic Tremor and Slip (ETS), involving transient deformations accompanied by emergent, low-frequency tremor occurs in subduction zones around the world. ETS events increase the shear stress on locked megathrusts and may potentially trigger damaging earthquakes. Despite the clear association of tremor and slip the physical relationship between them is unresolved. Tremor appears to result from slip on small asperities on the plate interface due to either creep on the surrounding fault, or stress increases ahead of the propagating slow-slip front. Previous studies of migrating slow slip events have not had sufficient spatial and temporal resolution to differentiate between these two models. To address this, we invert GPS data from the August 2009 ETS event in central Cascadia for the space-time evolution of fault slip-rate. We find a correlation in both space and time between tremor epicenters and the independently determined position of high fault slip-rate. This supports the first hypothesis that tremor asperities are loaded directly by slow slip, rather than by stress increases ahead of the slip front, and provides new insights into the mechanics of ETS.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011GL048714

Space-time correlation of slip and tremor during the 2009 Cascadia slow slip event

[1] Slow slip events (SSE) in many subduction zones incrementally stress the adjacent locked megathrust, suggesting that they could potentially either trigger or evolve into damaging earthquakes. We explore this with 2D quasi-dynamic simulations with rate-state friction, dilatancy, and coupled 1D pore-fluid and heat transport. Steady-state weakening friction allows transient slip to nucleate, but is inhibited by dilatant strengthening and destabilized by thermal pressurization. SSE spontaneously nucleate in Low Effective-Stress Velocity-Weakening (LESVW) regions. If the dimension of the LESVW is relatively small the SSE are trapped at its updip end, imparting a strong stress concentration in the locked zone. After several centuries SSE penetrate into the region of higher effective stress, where thermal pressurization eventually leads to dynamic rupture. For larger LESVW regions SSE tend to increase in length with time; ultimately higher slip speeds enhance thermal weakening, leading to dynamic instability within the SSE zone. In both cases the onset of the ultimate SSE is essentially indistinguishable from preceding events.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2012GL052811

Slow-slip evolves into megathrust earthquakes in 2D numerical simulations

Geophysical observations have shown that transient slow slip events, with average slip speeds v on the order of 10−8 to 10−7 m/s, occur in some subduction zones. These slip events occur on the same faults but at greater depth than large earthquakes (with slip speeds of order ∼ 1 m/s). We explore the hypothesis that whether slip is slow or fast depends on the competition between dilatancy, which decreases fault zone pore pressure p, and thermal pressurization, which increases p. Shear resistance to slip is assumed to follow an effective stress law τ = f (σ− p) ≡ f σ. We present two-dimensional quasi-dynamic simulations that include rate-state friction, dilatancy, and heat and pore fluid flow normal to the fault. We find that at lower background effective normal stress (σ), slow slip events occur spontaneously, whereas at higher σ, slip is inertially limited. At intermediate σ, dynamic events are followed by quiescent periods, and then long durations of repeating slow slip events. In these cases, accelerating slow events ultimately nucleate dynamic rupture. Zero-width shear zone approximations are adequate for slow slip events, but substantially overestimate the pore pressure and temperature changes during fast slip when dilatancy is included.

The Role of Thermal Pressurization and Dilatancy in Controlling the Rate of Fault Slip

Albany is a multiphysics code constructed by assembling a set of reusable, general components. It is an implicit, unstructured grid finite element code that hosts a set of advanced features that are readily combined within a single analysis run. Albany uses template-based generic programming methods to provide extensibility and flexibility; it employs a generic residual evaluation interface to support the easy addition and modification of physics. This interface is coupled to powerful automatic differentiation utilities that are used to implement efficient nonlinear solvers and preconditioners, and also to enable sensitivity analysis and embedded uncertainty quantification capabilities as part of the forward solve. The flexible application programming interfaces in Albany couple to two different adaptive mesh libraries; it internally employs generic integration machinery that supports tetrahedral, hexahedral, and hybrid meshes of user specified order. We present the overall design of Albany, and focus on the specifics of the integration of many of its advanced features. As Albany and the components that form it are openly available on the internet, it is our goal that the reader might find some of the design concepts useful in their own work. Albany results in a code that enables the rapid development of parallel, numerically efficient multiphysics software tools. In discussing the features and details of the integration of many of the components involved, we show the reader the wide variety of solution components that are available and what is possible when they are combined within a simulation capability.

/pdf/albany-using-component-based-design-to-develop-a-flexible-5b7cmyngri.pdf

Andrew M. Bradley

Papers

Dilatant strengthening as a mechanism for slow slip events

Space-time correlation of slip and tremor during the 2009 Cascadia slow slip event

Slow-slip evolves into megathrust earthquakes in 2D numerical simulations

The Role of Thermal Pressurization and Dilatancy in Controlling the Rate of Fault Slip

Albany: using component-based design to develop a flexible, generic multiphysics analysis code