The investigation into the physical mechanisms and behavior of geoelectric fields as potential precursors of earthquakes is a crucial yet challenging task in seismo-electromagnetics. Assuming the electrokinetic effect as a mechanism, we propose a novel numerical modeling algorithm based on the framework of the LAC GRTM to explore the behavior of geoelectric fields in stratified porous media during earthquake preparation. This algorithm incorporates two innovative aspects: (1) employing the Jordan decomposition to address the degeneracy problem encountered in quasi-stationary scenarios, and (2) utilizing digital linear filters for Hankel transforms to conduct the wavenumber integration. The accuracy of the algorithm in computing static displacements in an extremely low porosity half-space model is verified through comparison with analytical solutions. Numerical results from a half-space model demonstrate a strong consistency between the vertical electric field and filtration displacement. Notably, the maximum amplitude of the vertical electric field can reach approximately 2878.8 mV/km, which is detectable by receivers located at considerable distances from the epicenter. In addition, three distinct types of pre-earthquake electric fields are identified: (1) the localized electric field induced by slow P waves, (2) the interface electric responses, and (3) the direct converted electric field from the source. Results from a multi-layer porous model show the significant influence of the water table on the amplitude of the vertical electric field. This can be attributed to the influence of water saturation and permeability on the coefficient that combines the vertical electric field and the vertical filtration displacement. We also analyze the temporal evolution of different fields during earthquake preparation, which demonstrates a robust correlation between the temporal characteristics of pore pressure and that of the vertical electric field. This indicates that the vertical electric field offers an effective means to detect pore pressure evolution during earthquake preparation. Its sensitivity to the strike, rake, and dip angles of potential faults indicates its potential for determining fault geometry related to an impending earthquake. Furthermore, fluctuations in displacements enable a rough estimation of variations in stress and confining pressure. These observations underscore the advantage of integrating geoelectric precursor data with continuous displacement observations to enhance our understanding of earthquake preparation mechanisms.