Highlights: • A new method for LAD estimation from T-LiDAR in occluded areas is developed. • This method use kriging on LAD estimator (LAD-Kriging) and relies on TLS data only. • LAD-Kriging can be applied to all voxels whatever their reliability. • LAD-Kriging retrieves reliable estimates in not-explored and poorly-sampled voxels. • This method can be fitted to other unbiased estimators extracting metrics in point clouds. Abstract: Terrestrial Laser Scanning (TLS) has been used during the past decade to capture the complexity of 3D forest canopy structures, especially Leaf or Plant Area Density (LAD/PAD). TLS data, i.e. point cloud, can be divided into voxels to estimate the three-dimensional distribution of LAD/PAD. However, the combination effects of vegetation occlusion and shooting pattern of TLS scanners lead to a highly heterogeneous sampling, which limits the reliability of some local estimates, since several voxels are either not explored or poorly-explored by laser beams. In practice, recommendations vary regarding the minimum number of beams crossing voxels or the minimum path lengths required to provide reliable predictions. In addition, assigning a value to non-explored and poorly-explored voxels is still an open question. The present work proposes a new method, called LAD-kriging, to mitigate the impact of non-uniform sampling and to increase the accuracy of LAD estimates in non-explored and poorly-explored voxels. The method takes advantage of i) an unbiased LAD estimator of known variance, which was recently developed; ii) the spatial correlation of the LAD field, which derives from vegetation clumping. LAD-kriging computes kriging weights from mathematical derivations, which takes into account both spatial dependencies in the LAD field and the reliability of the estimate available in each voxel. It was evaluated through numerical experiments, which enabled to validate the algorithm and to evaluate its performance through comparison with true references. An example application to field data shows that such spatial correlations truly exist in the field and that LAD-kriging entails to reduce sampling errors (with respect to full resolution scanning). Although a real validation is impossible due to the lack of precise references at the voxel scale, this example allows to gain confidence for its application to field data. In our realistic numerical experiment, up to 25-30% of voxels can be explored by<100 beams, thereby leading to unreliable estimates of the LAD within these voxels and frequent large local errors. In cases where voxels were not explored at all, simple methods previously reported in earlier studies, such as ignoring occlusion or assigning "mean layer" values to unexplored voxels were inefficient (RMSE = 0.72 m(-1) and poor local accuracy). By contrast, LAD-kriging exhibited the smallest errors (RMSE = 0.48 m(-1)) and was locally more accurate. LAD-Kriging was also efficient for correcting voxels explored with few beams, i.e. unreliable, dropping RMSE from 0.92 to 0.42 m(-1) in these volumes. In practice, LAD-Kriging does not require defining a reliability threshold or any arbitrary parameterization, which is convenient for field applications and can be extended to other applications than LAD with TLS, such as UAS/VAS scanning, provided that unbiased estimators of known variance are available.