Hi @geolane93_KU,
Without seeing the data and having a better understanding of the purpose, it's difficult to give concrete recommendations. However, I do have a few thoughts that might help.
First, if you have ArcGIS Pro 3.0 or later, look into the Compare Geostatistical Layers tool. You can create various different EBK3D outputs and compare their cross validation statistics to see which are more accurate than others. Then can help choosing a subset size, transformations, and semivariogram models.
Second, a subset size of 20 sounds quite small to me, particular for the K-Bessel semivariogram. My experience is that you should use at least 50 points in each subset for a semivariogram model with so many parameters (and, usually, more than 100 is better).
Third, I would consider removing some of the surface points that may be playing too dominant of a role in the model. The problem is alleviated somewhat by using sectored neighborhoods, but the comparatively dense sampling at the surface is likely still negatively impacting subsurface predictions. In particular, I suspect that the estimated Elevation Inflation Factor (EIF) is being most affected here, and the EIF is an extremely important parameter for accurate results.
Fourth, if the jagged edges and artifacts are far away from the input points (like in the top or bottom corner of the 3D extent), then I would not worry too much about them. EBK (2D and 3D) often produces these kinds of artifacts when you extrapolate (predicting outside the input points), but it tends to be very stable when interpolating (predicting between the input points).