The main goal of this Ph.D. thesis was the development of an advanced numerical approach to simulate mass transport in microfluidic electrokinetic systems. Understanding of the electrokinetic aspects of mass transport in microfluidic systems has a great importance for a relatively new technique in separation science: capillary and microchip electrochromatography. This technique combines the advantages of high-performance liquid chromatography and capillary electrophoresis. Capillary electrochromatography uses the electroosmotic flow generated by a high, external voltage applied to drive the liquid phase through a porous medium. The general mathematical formulation of the problem of the electroosmotic flow includes the description of various phenomena of different nature: hydrodynamics, electrostatics, ion transport, adsorption and dissociation. In addition, actual microfluidic systems employed in electrochromatography demonstrate frequently the extremely intricate morphology associated with their porous structure impeding the numerical treatment. In this work, the electroosmotic flow is modelled by the iterative numerical solution of the coupled Poisson, Nernst-Planck, and Navier-Stokes equations. In order to realize a computational time required for large-scale simulations of mass transport in actual electrochromatographic systems, the developed numerical model was implemented at a parallel high-performance computer and then used to simulate various electrokinetic problems. Chapter 2 contains a brief general theoretical description of mass transport problems in polar liquids, which are typical for applications in electrochromatographic analysis. The behaviour of such liquids can drastically change when an external electric field is applied due to the presence of the electrical double layer at the solid-liquid interface. Chapter 3 describes mainly the lattice-Boltzmann formalism, an alternative approach in computational fluid dynamics, which allows easily to treat geometrically complex boundaries and which is inherently parallel. In this approach the fluid is modelled by particles moving on a regular lattice. At each time step the particles propagate to neighbouring lattice points and re-distribute their velocities in a local collision phase. This method is extended to electrohydrodynamic problems by incorporating in the model the Lorentz force arising from the interaction of electrical charges in the liquid with the applied electric field. In Chapter 4 the results of a number of simulations concerning various aspects of microfluidic electrokinetics are presented and discussed. They follow the description of the algorithm employed for the computer generation of confined random packings of spherical particles. That algorithm is based on an improved Jodrey-Tory procedure and allows to generate fixed beds of spheres with an arbitrary size distribution confined by an arbitrary container, as well as with periodic boundaries. The random sphere packings are further used as a model of particulate packed chromatographic columns. The presented numerical approach allows to obtain complete information concerning the spatial distribution in a modelled system of the flow velocity, electrical potential and species concentrations. In particular, the developed approach permits to evaluate the error related to the application of the apparent slip velocity boundary conditions to quantify differences between velocity fields obtained under different approximations concerning electrical boundary conditions, to study the effect of local variations in the chemical environment (caused by convection) on the surface charge density and final flow velocity field, to investigate the relation between the electroosmotic flow velocity and parameters of an electrokinetic system, such as the zeta-potential, solution concentration and applied electric field. In addition, the presented approach can be used to investigate the transient behaviour of simulated systems, such as the transient hydrodynamic dispersion in packed beds.