Abstract <p>In this work, we present a numerical model for the study of the generation of intense terahertz radiation in the interaction of two-color laser pulses of optical or infrared wavelength with argon, as well as the results of numerical modeling performed on its basis. The short terahertz pulses obtained in such a scheme can be used to generate strong slowly changing electric and magnetic fields of a given configuration. To simulate the macroscopic response of the plasma produced by ionization of a gas target by a bichromatic laser field a fully kinetic plasma model consisting of the Vlasov equations for the plasma distribution functions and Maxwell’s equations for the self-consistent electromagnetic field has been used. Ionization by the electric field is taken into account in the framework of the cascade mechanism using formulas for the tunnel ionization probability, which include the instantaneous value of the electric field strength. Our model allows for a self-consistent calculation of the current density of electrons generated during ionization, including the back reaction of coherent terahertz radiation on their post-ionization dynamics. The numerical code is based on the particle-in-cell method and employs state-of-the-art algorithms, such as Boris pusher to update the particle positions and velocities, finite difference time domain scheme for the electromagnetic fields, charge conservation method to satisfy the Gauss law for the electric field. Energy losses associated with ionization are accounted for by the introduction of an ionization current whose density is such that its scalar product with the electric field vector in each numerical cell is equal to the work of the electric field spent on ionization in that cell during one time step.</p>