Skip to main content
Geosciences LibreTexts

5.3 Seismic Refraction Jupyter Notebook

  • Page ID
    11982
  • Instructions

    • To use this page, read the text between each of the python code windows, then press RUN to execute the code in the box.
    • It is necessary to run each of the boxes in order.
    • Don't click restart in the cells unless you come back to the top of the page and start over.
    • If you modify the code, this modification will not remain after you leave or refresh this page.
    • You must run simulations/cells that require user input to completion. Make sure you see the Message informing you a cell is complete before continuing to the next cell. Here's why you must do this:
      • These cells may not be re-run until the simulation is completed, otherwise no output will be generated.
      • All following cells will not run until the user input cell is completed.
    • For Dropdown Menus: You do not need to rerun the cells with the menu to change your selection. However, you must rerun all of the following cells to implement the change.

    Be patient, some times it can take 1-2 minutes for the juypter kernel to start.

    An interactive example of single-layer seismic refraction.

    Key Questions: Consider these as you work your way through this page.

    • Holding the layer 1 thickness and velocity fixed, how does the angle of the incident wave (and refracted wave) change as you increase the layer 2 velocity?
    • How does the layer thickness affect the geometry of the ray path?
    • What happens if the velocity of the second layer is less that the first layer? Why does this happen?

    First, we import the necessary libraries.

    import numpy as np
    import matplotlib.pyplot as plt
    
    print('The Libraries are imported.')

    Now let's create our single-layer refraction model. This cell will prompt you to select 2 p-wave velocities and a layer thickness. It will then calculate the critical angle and plot the model. The model is created by plotting lines that reflect the desired geometry. Go ahead and run the cell, then consider the key questions.

    print('Typical p-wave velocities are between 1000 and 4000 m/s')
    print('Please enter the p-wave velocity of layer 1')
    v1 = int(input())
    print("Please enter the p-wave velocity of layer 2")
    v2 = int(input())
    print('Please enter the thickness of layer 1 between 0 and 500 m')
    h = int(input())
    while h <= 0 or h>500:
        print('Invalid thickness, enter a number between 0 and 500 m')
        h = float(input())
    
    h2 = 10 # m
    
    theta_ic = np.arcsin(v1 / v2)
    theta = np.pi/2 - theta_ic
    handle = display(None, display_id=True)
    bbox = dict(boxstyle ="round", fc = '1')
    while h > 0 and h <= 500:
        fig, ax = plt.subplots(figsize=(10,5))
    
        ax.set_xlabel('Distance from Source (m)')
        ax.set_ylabel('Depth (m)')
        ax.set_title('Seismic Refraction')
        ax.set_ylim(-h2,h+10)
        spacing = 10
        if h >= 200:
            spacing = 30
        ax.set_yticks(np.arange(-h2,h+spacing,spacing))
        ax.set_yticklabels(np.arange(h+h2,-spacing,-spacing))
        if v2 > v1:
            interface_distance = np.linspace(0,500,5)
            x1 = h * np.tan(theta_ic)
            xmax = interface_distance[-1]+2*x1 + 50
            ax.hlines(0,0,xmax,color='black')
            ax.hlines(h,0,xmax,color='black')
            ax.hlines(-h2,0,xmax,color='black')
            for i in range(len(interface_distance)):
                path = np.concatenate(((np.linspace(0,x1)*-np.tan(theta)+h),0*np.linspace(x1,interface_distance[i]+x1),(np.linspace(interface_distance[i]+x1,interface_distance[i]+(2*x1))-interface_distance[i]-x1)*np.tan(theta)))
                ax.plot(np.concatenate((np.linspace(0,x1),np.linspace(x1,interface_distance[i]+x1),np.linspace(interface_distance[i]+x1,interface_distance[i]+2*x1))),path,color='blue')
    
            ax.set_xlim(0,xmax)
            ax.text(xmax*1.05,h/2,'Layer 1 Velocity: '+str(v1)+' m/s',bbox=bbox)
            ax.text(xmax*1.05,-h2/2,'Layer 2 Velocity: '+str(v2)+' m/s',bbox=bbox)
            ax.text(xmax*1.05,h/4,'Critical Angle: '+'{:.0f}'.format(np.degrees(theta_ic))+' degrees',bbox=bbox)
        else:
            incoming_angle = np.radians(60)
            theta2 = np.arcsin((v2/ v1) * np.sin(incoming_angle))
            theta = np.pi/2 - theta2
            intercept = (h2-(h+h2))/(-np.tan(incoming_angle))
            xmax = (h+h2)/np.tan(incoming_angle) + 10
            ax.plot(np.linspace(0,intercept),np.linspace(0,intercept)*-np.tan(incoming_angle)+h,color='blue')
            ax.plot(np.linspace(intercept,xmax),(np.linspace(intercept,xmax)-intercept)*-np.tan(theta),color='blue')
            ax.hlines(0,0,xmax,color='black')
            ax.hlines(h,0,xmax,color='black')
            ax.hlines(-h2,0,xmax,color='black')
            ax.set_xlim(0,xmax)
            ax.text(xmax*1.05,h/2,'Layer 1 Velocity: '+str(v1)+' m/s',bbox=bbox)
            ax.text(xmax*1.05,-h2/2,'Layer 2 Velocity: '+str(v2)+' m/s',bbox=bbox)
        handle.update(fig)
        plt.close()
        print('To exit the simulation enter a number outside the range')
        print('Please enter the thickness of layer 1 between 0 and 500 m')
        h = int(input())
        
        
        
    print('The simulation is complete. Run the cell again to use different velocities.')

    Key Questions: Consider these as you work your way through this page.

    • Holding the layer 1 thickness and velocity fixed, how does the angle of the incident wave (and refracted wave) change as you increase the layer 2 velocity?
    • How does the layer thickness affect the geometry of the ray path?
    • What happens if the velocity of the second layer is less that the first layer? Why does this happen?