Nested numerical simulations of ionospheric plasma density structures associated with nonlinear evolution of the Rayleigh–Taylor (RT) instability in equatorial spread F (ESF) are presented. The numerical implementation of the nested model uses a spatial discretization with a C grid staggering configuration where normal velocities of ions and electrons are staggered one-half grid length from the density of charged particles. The advection of charged particles is computed with a fifth order accurate in space weighted essentially nonoscillatory (WENO) scheme. The continuity equation is integrated using a third-order Runge–Kutta (RK) time integration scheme. The equation for the electric potential is solved at each time step with a multigrid method. For the limited area and nested simulations, the lateral boundary conditions are treated via implicit relaxation applied in buffer zones where the density of charged particles for each nest is relaxed to that obtained from the parent domain. The high resolution in targeted regions offered by the nested model was able to resolve secondary RT instabilities, and to improve the resolution of the primary RT bubble compared to the coarser large domain model. The computational results are validated by conducting a large domain simulation where the resolution is increased everywhere.