Week of 2025.03.17

Here are my goals for this week:

• Work through math for stability conditions again and verify analytical results agree.
• Verify that Mario knows how to implement slip.
• Continue to noodle on the negative velocity feature. It appears to go away when time steps are small. Chat to Eric about a mitigation, or maybe it doesn't really matter.
• Review differential equation for slip. Calculate analytical solution for slip to compare to quail PDE solution.

1.0 Stability analysis reviewed

Note: The below calculations assume $L_{plug} >> s$.

Let's review the dynamics of our system:

$$ M \ddot{s} = A (p_0 + \Delta p(s)) - p_{atm} A - 2 \pi R L_{plug} \tau(s) $$

When the system is at equilibrium, the following is true:

$$ 0 = A (p_0 - p_{atm}) - 2 \pi R L_{plug} \tau_p $$

Equilibrium therefore implies

$$ p_0 = p_{atm} + \frac{2 L \tau_p}{R} $$

1.0.1 As a quick aside, lets review our steady state solution for our current parameters

$$ \begin{align} P_0 &= p_{atm} + \frac{2 L \tau_p}{R} \\ P_0 &= 1e5 [Pa] + \frac{2 * 100 [m] * 2e5 [Pa]}{10m} \\ P_0 &= 4.1e6 [Pa] \end{align} $$

Note, that we expect to see a linear pressure gradient of $p_{atm}$ -> $p_{0}$ across the plug, and then we expect to see the melt exist at the steady state pressure of $p_0$.

1.0.2 Stability conditions at steady state

To consider what values of $D_c$ and $s$ result in stability at equilibrium lets consider the dynamic equation at equilibirum. Enforcing the equlibrium conditions allows us to remove dependence on $p_0$ and $p_{atm}$. We get a new dynamics equation, at equilibrium of:

$$ M \ddot{s} = A \Delta p(s) + (\tau_p - \tau(s)) 2 \pi R L_{plug} $$

Intuitively, this means that if $ - A \Delta p(s) = (\tau_p - \tau(s)) 2 \pi R L_{plug}$ there will be no acceleration. At $s=0$ those two terms are equal ($ \Delta p(s) = 0$ and $\tau_s = \tau_p$).

The forces on both sides of the above equation are weakening forces. If (as s increases above 0), the $ A \Delta (p(s))$ grows more quickly is magnitude than $(\tau_p - \tau(s)) 2 \pi R L_{plug}$ that implies stability. Conversely, if the friction weakens more quickly than the pressure term, than the system will be unstable. Mathmatically, the system in unstable when:

$$ \begin{align} F_p(s) &= A \Delta p(s) \\ F_f(s) &= (\tau_p - \tau(s)) 2 \pi R L_{plug} \\ \frac{|\frac{d F_p}{ds}|}{|\frac{d F_f}{ds}|} &< 1 \end{align} $$

In English, this states that the system is unstable when the pressure force decreases more slowly than the frictional force decreases.

1.0.3 Fluid compressibility review

Let's review fluid compressibility and bulk modulus briefly before calculating the derivatives:

$$ \begin{align} K &= - V \frac{dP}{dV} \\ K &= - A * L \frac{dP}{dL * A} \\ K &= - L \frac{dP}{dL} \end{align} $$

1.0.4 Tau frictional force law (linear)

Recall that our linear model of $\tau$ is expressed like

$$ \tau(s) = \tau_p - (\tau_p - \tau_r) * \frac{s}{D_c} $$

1.0.5 Derivatives

$$ \begin{align} \frac{d F_P}{ds} &= A * \frac{dP}{ds} = \frac{-A * K}{L_{conduit}} \\ \frac{d F_f}{ds} &= 2 \pi R L_{plug} \frac{\tau_p - \tau_r}{D_c} \end{align} $$

1.0.6 $D_c$ as a function of $\tau_p$ and $\tau_r$ for instability

Instability conditions $$ \begin{align} \frac{|\frac{d F_p}{ds}|}{|\frac{d F_f}{ds}|} &< 1 \\ \\ \frac{\frac{A * K}{L_{conduit}}}{2 \pi R L_{plug} \frac{\tau_p - \tau_r}{D_c}} &< 1 \\ \\ \frac{A K D_c}{L_{conduit} 2 \pi R L_{plug} (\tau_p - \tau_r)} &< 1 \end{align} $$

So this results in an instability condition of

$$ \begin{align} D_c < \frac{L_{conduit} 2 \pi R L_{plug} (\tau_p - \tau_r)}{A K } \\ D_c < \frac{2 L_{conduit} L_{plug} (\tau_p - \tau_r)}{R K } \\ \end{align} $$

1.0.7 specific D_c for certain ICs test

If $\tau_p = 2e5$ and $\tau_r = 0$, $R = 10 [m]$, $L_{plug} = 100 [m]$, and $L_{conduit} = 900 [m]$ we get:

$$ D_c < \frac{2 * 900 [m] * 100 [m] * (2e5 [Pa])}{10 [m] * 1e9 [Pa] } \\ D_c < 3.6 [m] $$

This tells us that if D_c is less than 3.6 meters, the system will be unstable. If there is marginally more pressure than steady state, the plug will be ejected 3.6m before the loss in pressure will be greater than the reduction in friction.

1.0.8 Numerical test with $D_c < 3.6 [m]$

Let's set the parameters $D_c = 3.3 [m]$ and let's set the initial pressure in the metal equal to 4.5 MPa which is slightly higher than the 4.1 MPa of steady state.

from slip_imports import *

ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_v1")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

1.0.9 Numerical test with $D_c > 3.6 [m]$

Let's leave the parameters from above the same expect set $D_c = 3.7$. We now expect the system to quickly stabalize as the reduction in pressure happens more quickly than the fall off in pressure.

from slip_imports import *

ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_v4", p0=4.1)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

2.0.10 Changes to make analytical solution match quail solution

It turns out that the disagreement in analytical and numerical solutions were largely a result of me not considering the effect of viscous drag. By reducing viscous drag to $5e4$ I was able to largely solve this issue. The viscous drag term probably still causes the numerical model to be slightly more stable than my analytical solutions gives it credit for, but the analytical numerical comparison is much closer now.

2.0 Review of negative velocity

In the steady state solution, we see a negative velocity across the plug. This is unexpected given that we expect the frictional force to be equally weighted with the force from pressure. The negative velocity reduces the finer our timesteps indicating that this negative velocity is a function of numerical error.

In an email to me, Eric described the negative velocity as follows:

Regarding the negative velocity, I think it is simply the issue that I mentioned the other day, that the numerical steady state can be slightly different from the analytical steady state due to approximation of d/dx. The analytical steady state has a balance between the pressure gradient pushing in the +x direction and drag pushing in the -x direction (at least the way it's implemented now). If that pressure gradient is slightly different due to numerical approximation of dp/dx, then it's possible that the pressure gradient is a little bit too small and then we get a negative net force that pushes the plug down and creates a negative velocity.

What suprises me is that the numerical solution arrives at the same $P_0$ as our analytical solution, yet still introduces this negative velocity artifact. One of the challenges with the negative velocity artifact is that it reduces slip and therefore the system does not slip into the unstable configuration that we would otherwise expect.

2.0.1 Test Eric's hypothesis

As for negative velocity, try to test my hypothesis. From the output tau and p, calculate the pressure gradient term and the frictional drag term and check that drag force is a bit larger than pressure gradient force. If that's the case, then perhaps increase the pressure gradient a bit more in the initial condition?

Based on the below simulation, I believe Eric's hypothesis is false. I see that given the $\tau_p = 0.3 [MPa]$, the pressure gradient almost exactly matches our expectations. See above for the calculation expecting the pressure on the left side of the gradient to be 4.1 MPa.

For this simulation, $D_c = 100 [m]$ (large so that the system stays stable). The fact that $p_0$ lines up for the analytical and numerical simulation implies the pressure gradient is correct in the numerical sim.

from slip_imports import *

ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_v2", p0=4.1)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

2.0.2 Negative velocities at different number of time steps. Consider three examples of the steady solution calculuated at three different size of timesteps.

After chatting about this issue with Eric, we concluded that:

1. This negative velocity is likely a function of that absence on stick/slip logic. At some point we could consider implementing that.
2. I previously showed the negative velocity error decreased with decreases in mesh size. I am now deleting those graphs to reduce the size of this notebook.
3. Eric suggested I experiment turning viscous drag on inside the plug. The simulation below experiments with that. The conclusion from that experiment is that while adding viscous drag improves the situation (if we increase the viscosity substantially, it does not really solve the issue).

from slip_imports import *

ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=25, file_prefix="short_plug_v8", p0=4.1, max_velocity=0.1)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

3.0 Look at Mario's config file and see why it requires so many timesteps

The Courant Friedrich Lewy condition relates max wave speed with required time step.

$$\frac{u \Delta t}{\Delta x} \leq C_{max}$$

$C_{max}= 1$ because we are using an explicit time step method (according to wikipedia), so:

$$ \frac{\Delta x}{\Delta t} = u $$

In the below example, wave speed maxes out at 250 m/s and $\Delta x = 10m$, so the $\Delta t$ is on the order $\frac{1}{25}$.

Experimentally, I found that after fixing some weird behavior with the plug I was able to to reduce the number of timesteps to 2k for a 10 second simulation. That corresponds to 200 timesteps per second which is still below the max limit but a good deal faster than what we had been dealing with.

from slip_imports import *

ani = animate.animate_conduit_pressure("test_for_mario", iterations=50, file_prefix="tungurahua_rad_5_v14_conduit", viscosity_index=1, wall_friction_index=4, p0=4.1, x_min=-1000, max_speed_of_sound=1000, max_pressure=50, max_velocity=120)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

4.0 Solving the differential equation for slip

It is interesting to solve the differential equation for slip and compare to the results in our numerical simulation. This simulation ignores:

1. viscous drag
2. added mass effect

But it should give us a sanity check that the numbers we see in the above simulations (specifically 1.0.8 and 1.0.9 are legitimate).

From above, recall that our differential equation is (I am removing the simplifying assumption that $L_{plug} >> s$).

$$ \begin{align} M \ddot{s} &= A (p_0 + \Delta p(s)) - p_{atm} A - 2 \pi R L_{plug} \tau(s) \\ M \ddot{s} &= A * (p_0 - p_{atm} + \Delta p (s)) - \tau(s) * 2 \pi R (L_{plug} - s) \\ M \ddot{s} &= A * (p_0 - p_{atm} - \frac{K s}{L_{melt}}) - \tau (s) * 2 \pi R (L_{plug}-s) \\ \end{align} $$

Let's rewrite that as a two first order differential equations and then solve numerically for our parameters. Let's call slip $s_1$ and it's derivate $s_2$ $$ \begin{align} s_2 &= \dot{s_1} \\ \dot{s_2} &= \frac{A}{M} * (p_0 - p_{atm} - \frac{K s_1}{L_{melt}}) - \frac{\tau (s_1) * 2 \pi R (L_{plug}-s_1)}{M} \\ \end{align} $$

4.1 Updating the simulation to include added mass effect

To update the mass, we want to add the mass of the melt that begins moving to correctly slow down the acceleration. In practice the increase in the effective mass is a function of the ratio between $L_{melt}$ and $L_{plug}$. Specifically,

$$ M_{eff} = M_{plug} (1 + \frac{L_{plug}}{2 L_{melt}}) $$

In practice, the added mass effect appears to have minimal impact.

4.2 Updating lumped parameter model to include viscous drag

$$ \begin{align} F_{visc} &= \int_0^{L_m} 2 \pi R \tau_{visc} dz \\ F_{visc} &= \int_0^{L_m} 2 \pi R \frac{4 u \mu}{R} dz \\ F_{visc} &= 8 \pi \mu \int_0^{L_m} u dz \\ \end{align} $$

Let's assume that the velocity in the melt linearly increases from 0 at the bottom of the conduit to $V$ at the boundary of the melt and the plug $u = \frac{V * z}{L_{melt}}$. Integrating we get:

$$ F_{visc} = 4 \pi \mu L_{melt} V $$

In this case $V = \dot{s}$

So our updated differential equation becomes

$$ \begin{align} s_2 &= \dot{s_1} \\ \dot{s_2} &= \frac{A}{M_{eff}} * (p_0 - p_{atm} - \frac{K s_1}{L_{melt}}) - \frac{\tau (s_1) * 2 \pi R (L_{plug}-s_1)}{M_{eff}} - \frac{4 \pi \mu L_{melt} * s_2}{M_eff} \\ \end{align} $$

4.3 Updating lumped parameter model to include radiation damping

From Eric:

But the plot thickens! Mario and I worked through some other effects today and convinced ourselves that plug inertia is not important but another effect, known as radiation damping, is important. This is because the timescale of plug rupture and slip (estimated as ~1 s from Mario's analysis of the seismic data) is shorter than the timescale of waves propagating down to the bottom of the conduit and back (and those waves are what cause the pressure change in the magma to be spatially uniform, which was a key assumption in the lumped parameter model). Radiation damping is a solution that is valid in the short timescale limit, where as the plug slips, it sounds out an acoustic wave that propagates down through the magma. Solution to the wave equation provides a relation between the pressure decrease at the base of the plug and the slip velocity of the plug: Delta p = - rhocu, where rho = magma density, c = sound speed in magma. This is a ds/dt term in the lumped parameter model. you can compare that term to the d^s/dt^2 plug inertia term and show that the latter is smaller than the former, hence we can neglect plug inertia.

$$ \Delta p = - \rho * c *u$$

Some questions:

1. Is this radiation damping term already considered in the PDE model? I think the answer here is yes. The PDE model can models waves.
2. Given that the lumped parameter model is already has slower slip than the PDE model, will not adding the radiation damping term just slow it down even more? Yes, but perhaps as part of using a wave based model we will also decrease the drop in pressure based on the quasi static assumption.

4.3.1 Observations:

1. If I decrease the term for $\Delta p$ that is dependent on slip, I find that the steady state slip increases. Experimentally, reducing the the effect of slip on change in pressure by about 2 percent results in a numerical simulations that match.
2. In order to make the graphs match, I needed to flip the sign on the the radiation damping term and multiple by a constant of 0.5. So the numerical models only agree if the radiation damping results in higher pressure.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define constants (you can modify these values as needed)
K = 1e9      # Some constant (friction/damping)
L_melt = 900 # Melt length
tau_p = 2e5  # principle shear stress
tau_r = 0  # residual shear stress
R = 10      # Radius [m]
L_plug = 100 # Plug length

rho = 2.6e3  # Density

M_plug = R**2 * np.pi * L_plug * rho    # Mass
C = 600 # [m/s] Speed of sound in magma 

# Account for adding mass effect.
M_eff = M_plug * (1 + L_melt / (L_plug*2))
D_c = 3.3

# visocity
mu = 5e4

A = np.pi * R**2  # Cross-sectional area
p0 = 4.5e6  # Pressure at the top of the conduit
p_atm = 1e5  # Atmospheric pressure

def tau(s):
    #return tau_r - (tau_r - tau_p) * np.exp(s / D_c)

    if s < D_c:
        return tau_p - (tau_p - tau_r) * s / D_c
    else:
        return tau_r

# Define the system of differential equations
def system(state, t):
    s1, s2 = state  # s1 is position, s2 is velocity
    
    # ds1/dt = s2
    ds1_dt = s2
    
    # ds2/dt = (-A*K*s2)/(M*L_melt) + ((tau_p - tau_s)*2*pi*R*(L_plug - s1)*s1)/M
    ds2_dt =  (-A * K * s1) / (M_eff * L_melt) + \
            (A * (p0 - p_atm) / M_eff) + \
             - (tau(s1) * 2 * np.pi * R * (L_plug - s1)) / M_eff - \
             4 * np.pi * mu * L_melt * s2 / M_eff \
            #+ 0.5*(A * rho * C * s2 / M_eff) 
             
    
    return [ds1_dt, ds2_dt]

# Initial conditions
s1_0 = 0.0  # Initial position
s2_0 = 0.0  # Initial velocity
initial_state = [s1_0, s2_0]

# Time points
t = np.linspace(0, 45, 1000)  # Time from 0 to 10 with 1000 points

# Solve the differential equations
solution = odeint(system, initial_state, t)

# Extract solutions
s1 = solution[:, 0]  # Position
s2 = solution[:, 1]  # Velocity

# Plot results
plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
plt.plot(t, s1, 'b-', label='slip [m]')
plt.grid(True)
plt.xlabel('Time')
plt.ylabel('slip [m]')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(t, s2, 'r-', label='velocity [m/s]')
plt.grid(True)
plt.xlabel('Time')
plt.ylabel('velocity [m/s]')
plt.legend()

print(f"D_c: {D_c}")
plt.tight_layout()
plt.show()

# Optional: Print some key values
print(f"Slip: {s1[-1]:.4f}")
print(f"Velocity: {s2[-1]:.4f}")

D_c: 3.3

Slip: 3.9600
Velocity: -0.0000

Now, let's compare the numerical solution to our differential equation with the the results generated by quail. We will match the parameters up. Namely:

1. $P_0 = 4.5 Pa$ of pressure in the melt
2. $D_c = 3.3$
3. $\mu = 5e4$
4. $\tau_p = 2e5$
5. $\tau_r = 0$

folder = "slip_variable_addition"
file_prefix = "short_plug_v9"
N = 300
viscosity_index = 2
slip_time_series = []
times = []

for i in range(0, N):
    solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{i}.pkl")
    arhoA = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityA")]
    arhoWt = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityWt")]
    arhoWv = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityWv")]
    arhoC = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityC")]
    arhoM = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityM")]
    arhoF = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityFm")]

	# Get the value of the new state variable.
    rho_mix = arhoA + arhoWv + arhoM + arhoF
    momentum = solver.state_coeffs[:,:,solver.physics.get_momentum_slice()]
		
    # Get the value of the new state variable.
    slip_time_series.append((solver.state_coeffs[:,:,solver.physics.get_state_index("slip")]/rho_mix).ravel()[-1])
    times.append(solver.time)
# Plot results
plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
plt.plot(times, slip_time_series, 'b-', label='Quail PDE [m]')
plt.plot(t, s1, 'r--', label='Lumped parameter ODE[m]')
plt.xlabel("Time [s]")
plt.ylabel("Slip [m]")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

4.3 experimenting with dampening modes

For an ODE of the form:

$$ \begin{equation} m \ddot{s} = - k s - c \dot{s} \end{equation} $$

The dampening constant $\zeta = \frac{c}{2\sqrt{mk}}$

Rewriting out differential equation:

$$ \ddot{s} M_{eff} = A * (p_0 - p_{atm} - \frac{K s}{L_{melt}}) - \tau (s) * 2 \pi R (L_{plug}-s) - 4 \pi \mu L_{melt} * \dot{s} $$

Dropping the initial force from the pressure imbalence, assuming $s << L_{plug}$, and using the linear slip weakening model, we are left with:

$$ \ddot{s} M_{eff} = - (\frac{A K}{L_{melt}} - \frac{2 \pi R L_{plug} (\tau_p - \tau_r)}{D_c}) s - (4 \pi \mu L_{melt})\dot{s} $$

$$ \begin{align} \zeta &= \frac{c}{2\sqrt{mk}} \\ \zeta &= \frac{4 \pi \mu L_m}{2 \sqrt{M_{eff} * (\frac{AK}{L_m} - \frac{(\tau_p - \tau_r)2\pi R L_l}{D_c}})} \end{align} $$

When we calculate that below, we see that we are in the slightly underdamped regiment. I don't know how help this analysis is, because the linear aproximation I am using for the friction slip law is really only valid when $s < D_c$.

# Define constants (you can modify these values as needed)
K = 1e9      # Some constant (friction/damping)
L_melt = 900 # Melt length
tau_p = 2e5  # principle shear stress
tau_r = 0  # residual shear stress
R = 10      # Radius [m]
L_plug = 100 # Plug length

M_plug = R**2 * np.pi * L_plug * 2.6e3    # Mass

# Account for adding mass effect.
M_eff = M_plug * (1 + L_melt / (L_plug*2))
D_c = 3.3

# visocity
mu = 5e4

A = np.pi * R**2  # Cross-sectional area
p0 = 4.5e6  # Pressure at the top of the conduit
p_atm = 1e5  # Atmospheric pressure


def zi(D_c):
    return (4 * np.pi * mu * L_melt) / (2 * np.sqrt( M_eff * (A*K / L_melt - (tau_p - tau_r)*2*R*L_plug/D_c)))

zi(3.3)

np.float64(0.8837335780303277)