We present an efficient and robust numerical model for simulation of electrokinetic phenomena in porous networks over a wide range of applications including energy conversion, desalination, and 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 transport phenomena that can interact in complex and highly nonlinear ways in networks involving multiple pores with variable properties. Capturing these phenomena by direct simulation of the governing equations in multiple dimensions is prohibitively expensive. We present here a reduced order computational model that treats a network of many pores via solutions to 1D equations. Assuming that each pore in the network is long and thin, we derive a 1D model describing the transport in pore's longitudinal direction. We take into account the non-uniformity of potential and ion concentration profiles across the pore cross-section in the form of area-averaged coefficients in different flux terms representing fluid flow, electric current, and ion fluxes. Distinct advantages of the present framework include: a fully conservative discretization, fully bounded tabulated area-averaged coefficients without any singularity in the limit of infinitely thick electric double layers (EDLs), a flux discretization that exactly preserves equilibrium conditions, and extension to general network of pores with multiple intersections. By considering a hierarchy of canonical problems with increasing complexity, we demonstrate that the developed framework can capture a wide range of phenomena. Example demonstrations include, prediction of osmotic pressure built up in thin pores subject to concentration gradient, propagation of deionization shocks and induced recirculations for intersecting pores with varying properties.