For underground structures, a rough but reasonable simplification is pseudo-static deformation method. In this method, we apply seismic strain which can be calculated as the ratio of effective PGV (Peak Ground Velocity) to effective shear wave velocity.

$$ \gamma = \frac{PGV_e}{VS_e}  $$


Effective PGV can be multiplied with depth dependent reduction factors (see FHWA-NHI-10-034) and maximum shear wave velocity obtained from geophysical tests with almost zero strain can be converted to effective shear wave velocity based on recommendations of FHWA or Eurocode 8.
A simple Python code can be written to implement lateral deformation profile in Plaxis to simulate seismic loading. If you locate this Python file inside the Bentley folder (< PLAXIS installation folder >\pytools\input) this can be directly called from Plaxis Input.

"""
Seismic Deformations for Plaxis by Berk Demir / https://github.com/berkdemir
Locate this file in <PLAXIS installation folder>\pytools\input and call it from Plaxis - Expert - Run Python tool
"""

import easygui
from plxscripting.easy import *
# If script is used OUTSIDE of Plaxis.
"""
localhostport_i = 10000 # Local host port id can be different.
password = "password" # Your password should be here.
s_i, g_i = new_server("localhost", localhostport_i, password=password)
"""
# If we use the script as a Plaxis tool.

s_i, g_i = new_server()

def get_borehole_layers(borehole):
    """ reads the borehole information to collect soillayer thickness
        information and returns a dictionary per layer top-down """
    borehole_layers = []
    for soillayer in g_i.Soillayers:
        for zone in soillayer.Zones:
            if (zone.Borehole.value) == borehole:
                borehole_layers.append({"name": soillayer.Name.value,
                                        "top": zone.Top.value,
                                        "bottom": zone.Bottom.value,
                                        "thickness": zone.Thickness.value
                                        }
                                       )
    return borehole_layers

def get_xmin_xmax():
    """ gets the xmax and xmin of the model."""
    point_list = []
    
    for i in g_i.SoilContour:
        point_list.append(i)
    
    xmin = point_list[0].x.value
    xmax = point_list[1].x.value
    return xmin, xmax

bh = g_i.Boreholes[0]

borehole_layers = get_borehole_layers(bh)

top = borehole_layers[0]["top"]
bottom = borehole_layers[-1]["bottom"]
depth = top-bottom

title = "Seismic Deformation Application by Berk Demir"
msg = "Please enter seismic parameters."
fieldNames = [
    "Peak Ground Velocity (cm/sec)",
    "Reduction in PGV due to depth",
    "Maximum Shear Wave Velocity (m/s)",
    "Ratio of Effective to Maximum Shear Wave Velocity"
    ]

fieldValues = easygui.multenterbox(msg, title, fieldNames)

def_choice = ["Triangular","Z-Shape"]

select_def_type = easygui.buttonbox("Select deformation type", "Seismic Deformation Type", def_choice)

PGV, PGV_Red, VS, VS_Red = [float(item) for item in fieldValues]

xmin, xmax = get_xmin_xmax()

strain = PGV * PGV_Red * 0.01 / (VS*VS_Red)

if select_def_type == "Triangular":
    deformation = strain * depth
    LD_Left = g_i.linedispl((xmin,top),(xmin,bottom))[-1]
    LD_Right = g_i.linedispl((xmax,top),(xmax,bottom))[-1]
    LD_Top = g_i.linedispl((xmin,top),(xmax,top))[-1]
    LD_Left.setproperties("Displacement_x","Prescribed","Displacement_y","Free","Distribution","Linear","ux_start",deformation,"ux_end",0)
    LD_Right.setproperties("Displacement_x","Prescribed","Displacement_y","Free","Distribution","Linear","ux_start",deformation,"ux_end",0)
    LD_Top.setproperties("Displacement_x","Prescribed","Displacement_y","Fixed","Distribution","Uniform","ux_start",deformation)
    
elif select_def_type == "Z-Shape":
    deformation = strain * depth * 0.5
    LD_Left = g_i.linedispl((xmin,top),(xmin,bottom))[-1]
    LD_Right = g_i.linedispl((xmax,top),(xmax,bottom))[-1]
    LD_Top = g_i.linedispl((xmin,top),(xmax,top))[-1]
    LD_Bottom = g_i.linedispl((xmin,bottom),(xmax,bottom))[-1]
    LD_Left.setproperties("Displacement_x","Prescribed","Displacement_y","Free","Distribution","Linear","ux_start",deformation,"ux_end",-deformation)
    LD_Right.setproperties("Displacement_x","Prescribed","Displacement_y","Free","Distribution","Linear","ux_start",deformation,"ux_end",-deformation)
    LD_Top.setproperties("Displacement_x","Prescribed","Displacement_y","Fixed","Distribution","Uniform","ux_start",deformation)
    LD_Bottom.setproperties("Displacement_x","Prescribed","Displacement_y","Fixed","Distribution","Uniform","ux_start",-deformation)

To determine the effective shear wave velocity, I have developed a simple methodology and will present this as a part of a paper in the near future (hopefully in WTC 2022). I will not go into detail too much, but you can understand from the code that it is a iterative process.

  • Calculate seismic shear strain using reduced PGV based on the depth of tunnel and maximum shear wave velocity.
  • Ignoring any static strain, calculate the G/Gmax based on either one of the presented approaches.
  • Calculate shear modulus reduction ratio as square root of the G/Gmax.
  • Using the calculated ratio, calculate effective shear wave velocity.
  • Calculate new seismic shear strain using reduced PGV and new effective shear wave velocity.
  • Continue until reasonable difference is obtained between Vsi+1 and Vsi.

Here is the simple code that calculates the effective shear wave velocity:

def EFS(Ground_Type, PGV, PGV_Red, VS, PI=20, OCR=1, Eff_Pressure=300):
    """
    Effective Shear Wave Velocity Using Darendeli (2001) and Schnabel (1973) by Berk Demir / https://github.com/berkdemir
    
    Inputs:
        Ground_Type: Either "Soil" or "Rock"
        PGV: Peak Ground Velocity (cm/sec) 
        PGV_Red: PGV reduction ratio with depth
        VS: Maximum shear wave velocity (m/s)
        PI: Plasticity index in percent. (For Soil only.)
        OCR: Overconsolidation Ratio. (For Soil only.)
        Eff_Pressure: Effective pressure (kPa) (For Soil only)
        
    Returns:
        VS_Red: Reduction ratio for maximum shear wave velocity.
    
    Notes:
        Also prints a sentence with reduction ratio and effective shear wave velocity.
        
    Example:
        For Soils: Vs_Red = EFS(Ground_Type = "Soil", PGV = 70, PGV_Red = 0.8, VS = 300, PI=20, OCR=1, Eff_Pressure=300)
        For Rocks: Vs_Red = EFS(Ground_Type = "Rock", PGV = 70, PGV_Red = 0.8, VS = 800)
    """
    
    PGV_eff = PGV * PGV_Red * 0.01  # m/s
    Shear_Modulus_Reaction = 0.7  # initial value
    VS_Red = pow(Shear_Modulus_Reaction, 0.5)

    if Ground_Type == "Soil":
        Strain_Ref = (0.0352 + 0.001 * PI * pow(OCR, 0.3246)) * \
            pow(Eff_Pressure / 100, 0.3483) / 100
        VS_Red = pow(Shear_Modulus_Reaction, 0.5)
        for i in range(0, 20):
            seismic_shear_strain = PGV_eff / (VS_Red * VS)
            Shear_Modulus_Reaction = 1 / \
                (1 + pow(seismic_shear_strain / Strain_Ref, 0.919))
            VS_Red = pow(Shear_Modulus_Reaction, 0.5)
        print("Reduction ratio of shear wave velocity is", round(VS_Red, 2), "and effective shear wave velocity is", round(
            VS*VS_Red, 0), "m/s using Darandeli (2001) G/Gmax curves for soils.")

    elif Ground_Type == "Rock":
        for i in range(0, 20):
            seismic_shear_strain = PGV_eff / (VS_Red * VS)
            Shear_Modulus_Reaction = -3.642 * \
                pow(seismic_shear_strain, 0.02553) + 3.784
            VS_Red = pow(Shear_Modulus_Reaction, 0.5)
        print("Reduction ratio of shear wave velocity is", round(VS_Red, 2), "and effective shear wave velocity is", round(
            VS*VS_Red, 0), "m/s using Schnabel (1973) G/Gmax curve for rocks.")

    return VS_Red

if __name__ == "__main__":
    VS_Red = EFS(Ground_Type, PGV, PGV_Red, VS, PI, OCR, Eff_Pressure)
    VS_eff = VS * VS_Red