Univariate cubic L1 spline fits have been successful to preserve the shapes of 2D data with abrupt changes. The reason is that the minimization of L1 norm of the data is considered, as opposite to L2 norm. While univariate L1 spline fits for 2D data are discussed by many, bivariate L1 spline fits for 3D data are yet to be fully explored. This thesis aims to develop bicubic L1 spline fits for 3D data approximation. This can be achieved by solving a bi-level optimization problem. One level is bivariate cubic spline interpolation and the other level is L1 error minimization. In the first level, a bicubic interpolated spline surface will be constructed on a rectangular grid with necessary first and second order derivative values estimated by using a 5-point window algorithm for univariate L 1 interpolation. In the second level, the absolute error (i.e. L1 norm) will be minimized using an iterative gradient search. This study may be extended to higher dimensional cubic L 1 spline fits research.