Abstract: Airborne Light Detection and Ranging (LiDAR) technology has widely been used to efficiently acquire terrain data over large areas, particularly providing data sources for the generation of high-resolution Digital Elevation Models (DEMs). However, little was known about the errors in the interpolation of airborne LiDAR point clouds for topographically complex areas, thereby resulting in the less application of airborne LiDAR in the earth surface process. In this study, the errors of six commonly-used DEM interpolations were assessed using the airborne LiDAR point clouds acquired from a topographically complex area in the gullied Loess Plateau, China. Six DEM interpolations included the Inverse Distance Weighted (IDW), Kriging, Spline, Natural Neighbor (NN), Triangulated Irregular Network (TIN), and Trend. Firstly, four parameters were optimized, including IDW, Kriging, Spline, and Trend. Secondly, the optimized algorithms were applied to produce DEMs. Lastly, the errors of DEMs were quantitatively evaluated using mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R2). Results showed that: 1) The optimal values of IDW weights and searching points were 1 and 12, respectively. The optimal parameters for Kriging included non-orientation, Gauss function, and 12 searching points. Spline performed best in regular spline functions and 32 of searching points. Nevertheless, Trend was unsuitable for the topographically complex area, due to the produced DEMs with meter-level errors. 2) In terms of quality, IDW, Kriging, NN, and TIN produced relatively sound DEMs at the 1%-25% subsampling rate (1-19 points /m2), whereas, unreasonable outliers were found in the DEMs produced by Spline. Additionally, there were similar spatial patterns in the DEMs produced by IDW, Kriging, NN, TIN, and Spline, when the subsampling rate was 50%-100% (39-77 points /m2). 3) Excellent relations (R2>0.99) were found between the elevation measurements and the DEMs produced using point clouds of different subsampling rates. The average MAE and RMSE of produced DEMs firstly decreased rapidly, and then stabilized, as the point density increased, demonstrating that the reduction of interpolation errors varied slowly, as the point density reached a certain level. TIN produced the lowest error at a 1% subsampling rate (1 points/m2), with the MAE and RMSE of 0.208 and 0.298 m, respectively. At the 5%-12.5% subsampling rate (4-10 points/m2), NN produced the lowest error, where the MAE and RMSE were 0.170-0.175 m and 0.259-0.262 m, respectively. At the >25%-100% subsampling rate (>19-77 points/m2), Spline yielded the lowest error with the MAE and RMSE of 0.164-0.168 and 0.249-0.255 m, respectively. More importantly, the interpolation errors for steep areas were considerably higher than those for gently-sloping areas. The errors for steep areas decreased markedly, while those for gently-sloping areas changed slightly, as the point density increased. NN and TIN were the most suitable interpolation for gently- and steep-sloping areas at the point density of 1-19 points/m2, with the MAE of 0.049-0.171 and 0.062-0.776 m, respectively. Spline yielded the lowest interpolation errors for both steep- and gently-sloping areas with the MAE of 0.010-0.123 and 0.051-0.593 m, respectively when the point density was between 39-77 points/m2. The findings can provide promising potential support to the earth surface process, thereby generating high-resolution DEMs for the topographically complex areas.