AbstractThis paper introduces a method for the optimal generation of multivariate intensity measure (IM) maps representing various aspects of the intensity of an earthquake over a region that can be used as input for regional hazard and loss-estimation analyses. The proposed method is an extension of the single-variate hazard quantization (HQ) methodology for the effective sampling of IM maps with a single intensity measure. The use of multiple intensity measures to describe an earthquake enables a more comprehensive representation of the event, with richer information. For instance, to conduct accurate regional seismic loss analyses of a portfolio of structures with different fundamental periods, maps of spectral accelerations at different periods (multivariate IM maps) are necessary. With this in view, the multivariate HQ method is proposed for the effective generation of a relatively small set of multivariate IM maps that can capture the seismicity of a region. A case study of the seismic region of Charleston, South Carolina, is presented. The sample space is represented by a large set of pairs of IM maps representing the spectral acceleration at two different periods, with each pair representing a specific earthquake event. The simulated pairs of IM maps carry the stochastic characteristics of the associated seismic events, such as the spatial correlation of the ground motions and the cross correlation between different intensities at different sites in the region. The ability of the proposed method to provide an accurate estimation of the hazard curve for both intensity measures and to correctly capture the spatial autocorrelations and the spatial cross correlation among the two intensity measures was investigated quantitatively. Experiments were conducted to demonstrate the robustness of the method and the effect of the sample size on the performance of the method.