A postprocessor has been developed to calculate space/time distributions of electrokinetic potentials resulting from histories of underground conditions (pressure, temperature, vapor saturation, concentrations of dissolved species, flow rate, etc.) computed by unsteady multidimensional geothermal reservoir simulations. Electrokinetic coupling coefficients are computed by the postprocessor using formulations based on experimental work reported by Ishido and Mizutani <1981>. The postprocessor was applied to both numerical modeling of natural self-potential (SP) anomalies in geothermal fields and production-induced SP changes. The essential features of the SP anomalies and SP changes which have actually been observed in real geothermal fields are reproduced reasonably well in these calculations for both single-phase (liquid) and two-phase (vapor/liquid) geothermal reservoirs. The postprocessor was also applied to predicting the magnitude of electrical earthquake precursory signals caused by dilatant strains taking place within a narrow vertical fault zone. If the rate of porosity increase is high enough to induce substantial pore pressure decrease, two-phase flow develops in the fault zone and produces observable signals at the ground surface if a high conductivity channel connects the near-surface region and the deep fault zone. ¿ 1999 American Geophysical Union