It has long been recognized that the mean provides an inadequate summa
ry whereas the set of quantiles can supply a more complete description
of a sample. We introduce bivariate quantile smoothing splines, which
belong to the space of bilinear tensor product splines, as nonparamet
ric estimators for the conditional quantile functions in a two-dimensi
onal design space. The estimators can be computed by using standard li
near programming techniques and can further be used as building-blocks
for conditional quantile estimations in higher dimensions. For modera
tely large data sets, we recommend penalized bivariate B-splines as ap
proximate solutions. We use real and simulated data to illustrate the
methodology proposed.