The one-dimensional fast Fourier transform (FFT) has been applied extensively to simulate Gaussian random wave elevations and water particle kinematics. The actual sea elevations/kinematics exhibit non-Gaussian characteristics that can be represented mathematically by a second-order random wave theory. The elevations/kinematics formulations contain frequency sum and difference terms that usually lead to expensive time-domain dynamic analyses of offshore structural responses. This study aims at a direct and efficient two-dimensional FFT algorithm for simulating the frequency sum terms. For the frequency-difference terms, inverse FFT and forward FFT are implemented, respectively, across the two dimensions of the wave interaction matrix. Given specified wave conditions, the statistics of simulated elevations/kinematics compare well with not only the empirical fits but also the analytical solutions based on a modified eigenvalue/eigenvector approach, while the computational effort of simulation is very economical. In addition, the stochastic analyses in both time domain and frequency domain show that, attributable to the second-order nonlinear wave effects, the near-surface Morison force and induced linear oscillator response are more non-Gaussian than those subjected to Gaussian random waves.