A biobjective PDE constrained shape optimization problem where the reliability and the cost of a ceramic component is optimized simultaneously and for which gradient information is available is considered. The probability of failure of the component under tensile load is modeled via a probabilistic Weibull-type model, while it is assumed that the cost is proportional to the volume of the component. Three different classes of gradient-based biobjective optimization methods are suggested and compared at 2D test cases. For the numerical implementations a first discretize then optimize strategy is utilized. The implementations benefit from efficient gradient computations using adjoint equations. The first class of gradient-based optimization methods that are considered are biobjective decent methods, where a biobjective steepest descent algorithm and the weighted sum method are investigated. The results of these methods yield approximations of the Pareto fronts that nicely display the trade-off between reliability and cost and give rise to innovative shapes that compromise between these conflicting objectives. For the second class of methods a novel approach for the efficient and reliable approximation of the Pareto front of sufficiently smooth unconstrained biobjective optimization problems is suggested and compared with a predictor-corrector method from the literature. For the novel method optimality conditions of the weighted sum scalarizations are utilized to describe (parts of) the Pareto front as a parametric curve, parameterized by the scalarization parameter (i.e., the weight in the weighted sum scalarization). Its sensitivity w.r.t. parameter variations are described by an ordinary differential equation (ODE). Starting from an arbitrary initial Pareto optimal solution, e.g., a solution of the first class of methods, the Pareto front can then be traced by numerical integration. An error analysis based on Lipschitz properties is provided and an explicit Runge-Kutta method for the numerical solution of the ODE is suggested. This method is then validated on biobjective test functions and applied on the given biobjective shape optimization problem. The third class of gradient-based biobjective optimization methods are surrogate model based methods. In particular, „Efficient Global Optimization“ (EGO) is considered which utilizes a surrogate model that is constructed via Gaussian random fields and is a global optimization method. The Gaussian model is also called Kriging model. In this work, we also consider Gradient Enhanced Kriging (GEK) models where gradient information is incorporated in the Kriging model. We choose EGO with Kriging and GEK to compute benchmark solutions for the given biobjective shape optimization problem since EGO is commonly used in applications. These benchmark results are then compared with the results of the other classes of methods.