We present an efficient and robust numerical model for the simulation of electrokinetic phenomena in porous media and microstructure networks considering a wide range of applications including energy conversion, deionization, and microfluidic-based lab-on-a-chip systems. Coupling between fluid flow and ion transport in these networks is governed by the Poisson Nernst Planck Stokes equations., These equations describe a wide range of phenomena that can interact in a complex fashion when coupled networks involving multiple pores with variable properties. Capturing these phenomena by direct simulation of the governing equations in multidimensions is prohibitively expensive. We present here a reduced order model that treats a network of many pores via solutions to 10 equations. Assuming that each pore in the network is long and thin, we derive a ID model describing the transport in the pore's longitudinal direction. We take into account the cross-sectional nonuniformity of potentiaL and ion concentration fields in the form of area averaged coefficients in different flux terms representing fluid flow, electric current, and ion fluxes. These coefficients are obtained from the solutions to the Poisson Boltzmann equation and are tabulated against dimensionless surface charge and dimensionless thickness of the electric double layer (EDL). Although similar models have been attempted in the past, distinct advantages of the present framework include a fully conservative discretization with zero numerical leakage, fully bounded area-averaged coefficients without any singularity in the limit of infinitely thick EDLs, a flux discretization that exactly preserves equilibrium conditions, and extension to a general network of pores with multiple intersections. In part II of this two-article series, we present a numerical implementation of this model and demonstrate its applications in predicting a wide range of electrokinetic phenomena in microstructures.