AbstractOn Mars, dust devils play an important role in injecting dust grains into the atmosphere. The exact amount that they contribute to the dust budget of the atmosphere is yet not clearly known. In this study, we model the spatial distribution of dust concentration within a steady state Martian dust devil for the first time. We numerically solve the equations of motion for dust particles to determine their velocity inside a dust devil, 10 m wide and 1,000 m tall, and consequently determine the dust loading using the continuity equation. We consider an initial wind profile, which is dependent on the circulation strength of the vortex (Γ) and viscosity of the air (ν). Our simulations indicate a maximum concentration of ∼1,400 cm−3 near the surface and at the boundary of the vortex. The larger size particles are lifted to lower heights. The radial and tangential particle velocities peak at ∼60 m, while the vertical velocity peaks at ∼100 m. A higher circulation strength (Γ), leads to a higher loading of dust, whereas a change in the air viscosity (ν) does not have a significant effect on the dust loading inside the steady state dust devil. From the simulated dust distribution in our vortex, the estimations of a dust flux of ∼5 × 10−5 kgm−2s−1, a total optical depth of 0.2 and a near‐surface heating rate of , are consistent with observations. Our calculations can provide useful inputs to study the effect of dust devils on boundary layer processes.