def _calc_Gmax(Vs, rho):
    """
    Calculate initial stiffness of each soil layer.

    Parameters
    ----------
    Vs : numpy.ndarray
        1D array of Vs of layers. Unit: m/s.
    rho : numpy.ndarray
        1D array of mass density of layers. Unit: kg/m^3.

    Returns
    -------
    Gmax : numpy.ndarray
        1D array of initial stiffness. Unit: Pa
    """
    Gmax = rho * Vs ** 2
    return Gmax
