Organic electrochemical transistors (OECTs) are of significant research interest as building blocks for the realization of synaptic behavior, bioelectronics, and neuromorphic applications. This work presents a numerical method considering the Poisson-Boltzmann statistics to reproduce associated charge densities and distributions depending on pH values as well as the electrolyte behavior, thus enabling numerical modeling and TCAD simulation of OECTs. We have fabricated and characterized OECTs based on PEDOT:PSS as a channel material. The proposed model does not take non-faradaic processes into account, but yet it reveals important properties of the device working mechanism and shows a good agreement with the measured data of OECTs.