We use recent results on statistical analysis of seismicity to present a robust method for comprehensive detection and analysis of earthquake clusters. The method is based on nearest‐neighbor distances of events in space‐time‐energy domain. The method is applied to a 1981–2011 relocated seismicity catalog of southern California having 111,981 events with magnitudes m ≥ 2 and corresponding synthetic catalogs produced by the Epidemic Type Aftershock Sequence (ETAS) model. Analysis of the ETAS model demonstrates that the cluster detection results are accurate and stable with respect to (1) three numerical parameters of the method, (2) variations of the minimal reported magnitude, (3) catalog incompleteness, and (4) location errors. Application of the method to the observed catalog separates the 111,981 examined earthquakes into 41,393 statistically significant clusters comprised of foreshocks, mainshocks, and aftershocks. The results reproduce the essential known statistical properties of earthquake clusters, which provide overall support for the proposed technique. In addition, systematic analysis with our method allows us to detect several new features of seismicity that include (1) existence of a significant population of single‐event clusters, (2) existence of foreshock activity in natural seismicity that exceeds expectation based on the ETAS model, and (3) dependence of all cluster properties, except area, on the magnitude difference of events from mainshocks but not on their absolute values. The classification of detected clusters into several major types, generally corresponding to singles, burst‐like and swarm‐like sequences, and correlations between different cluster types and geographic locations is addressed in a companion paper. Key Points Earthquake clusters are identified in southern California Accuracy and stability of detection is tested using ETAS model Several new cluster features are reported