Pore network models are efficient tools for upscaling flow and transport properties in porous media. This work introduces a new formal mathematical derivation of the discrete equations governing solute transport in a pore network model. A double Laplace transform technique is applied by enforcing mass flux continuity along with the interfaces between pores and throats. A non-local semi-analytical formulation results. Solutions are given as a sum of convolution products with time-dependent and exponentially decaying local Péclet-dependent infinite series. Continuous concentration profiles along throats are calculated analytically, a posteriori, from time-dependent numerically simulated concentrations in pores. The upwind and central-difference schemes of the widely used mixed-cell method are found to be equivalent to the asymptotic form of this new formulation for the advective and diffusive dominant regimes, respectively. Therefore, the validity range of these static methods is established. The model was compared to the delay differential equations approach, a newly derived analytical solution, and mixed-cell methods on idealized one-dimensional networks eliminating topological disorder. Concentrations in pores are best reproduced when transport in the throats is not neglected unlike for the mixed-cell method leading to early breakthrough and first-order moment. An efficient numerical scheme truncating the encoded memory effects in the convolution kernels is introduced. This paves the way for the model application to realistic networks extracted from digital rock images. We caution against using static formulations as the error can be very large locally before attending a steady-state.