Geostatistical conditional simulations are used to get a family of non-smoothed three-dimensional (3D) images of the St. Lawrence krill aggregation from echointegration data at 38 and 120 kHz on a systematic grid of transects. These maps respect the inherent patchiness of krill. Their ensemble gives an histogram of the krill density estimates at any point of the 3D grid, not only the mean density as in kriging. The spatially consistent simulations are conditional to match the observations and their histogram and variogram. The 3D problem is by-passed by transposing to a short series of 2D ones, which simplifies the modelling of the autocorrelation function and has the additional advantage of reducing by a factor of 5 the number of points to simulate. The spatially structured krill density profile is summarised by the first few factors of a principal component analysis, where the variables are the volume backscattering strength at 120 kHz for the integrated vertical bins along the profiles and the observations are the different profiles. The principal component scores are simulated over a 2D-grid and then used to reconstruct the full 3D-image. The method adequately reproduces the histogram, variogram and mean profile of krill density, and cross-validations replicate the observations. It is as precise as an alternative kriging approach to estimate the mean density, but has the additional advantages of requiring less computing time for our particular application. These conditional simulations enable estimating the probability density function of the krill density at any point, an essential information for predator-prey interactions and other ecosystem studies.