This contribution presents a method for the numerical determination of the steady-state response of complex charged porous media to pressure, salt concentration and electric potential gradients. The Pore Network Model (PNM), describing the porosity as a network of pores connected by channels, is extended to capture electrokinetic couplings which arise at charged solid-liquid interfaces. This allows us to compute the macroscopic fluxes of solvent, salt and charge across a numerical sample submitted to macroscopic gradients. On the channel scale, the microscopic transport coefficients are obtained by solving analytically (in simple cases) or numerically the Poisson-Nernst-Planck and Stokes equations. The PNM approach then allows us to upscale these transport properties to the sample scale, accounting for the complex pore structure of the material via the distribution of channel diameters. The Onsager relations between macroscopic transport coefficients are preserved, as expected. However, electrokinetic couplings combined with the sample heterogeneity result for some macroscopic transport coefficients (e.g. permeability or electro-osmotic coefficient) in qualitative differences with respect to their microscopic counterparts. This underlines the care that should be taken when accounting for transport properties based on a single channel of average diameter.