The calculation of earthquake hypocenters requires careful treatment, particularly when prior knowledge of the study area is limited. The prior knowledge, such as wave velocity and data noise, is often assumed to be known in earthquake location algorithms. Such assumptions can greatly simplify the inverse problem but are less general than nonlinear approaches. A nonlinear treatment is of particular importance when the uncertainty quantification of locations is of interest. We present a nonlinear multiple‐earthquake location method that is applicable when little prior knowledge of the area exists. Efficient Markov chain Monte Carlo (MCMC) sampling is employed in conjunction with a hierarchical Bayesian model that treats earthquake hypocenter parameters, as well as P‐wave velocity, ratio in P‐/S‐wave velocities, and P‐ and S‐data noise standard deviations as unknown. Hypocenters for multiple earthquakes are located concurrently to provide sufficient constraints for the parameter’s P‐wave velocity, ratio in P‐/S‐wave velocity, and P‐ and S‐data noise standard deviations, which are shared among events. The algorithm is applied to simulated and field data. With field data, 47 event hypocenters are located in 1 yr of data from 10 sensors in the Canadian Rocky Mountain trench. To analyze the probabilistic solutions, we compare single‐earthquake and multiple‐earthquake locations for the 47 events and find that the multiple‐earthquake location produces better‐constrained solutions when compared with the single‐event case. In particular, depth uncertainties are significantly reduced for the multiple‐earthquake location. The algorithm is inexpensive, considering that it is based on an MCMC approach and highly objective, requiring little practitioner choice for tuning.