A computer implemented method of generating an unstructured grid of a model of a subsurface reservoir based on input data defining internal boundary geometry descriptions of the reservoir model, near internal boundary grid cell points being generated from the intersections of suitable internal boundary intersected circles of adjacent boundary vertex points, an internal boundary being indicated by a discontinuity within the reservoir model, the method comprising the computer implemented steps of: (a) refining the internal boundary geometry data to produce a partial set of geometry data points to represent the internal boundary geometry taking into account a near internal boundary grid spacing requirement, said refining comprises selecting a partial set of geometry data points where the distance between adjacent geometry data points satisfies the near internal boundary grid spacing requirements; (b) building lines by connecting adjacent geometry data points from the step above to form internal boundary line segments ; (c) constructing internal boundary intersected circles at each of the geometry data points, wherein the radius of each intersected circle is selected so that it is greater than a half-length of a corresponding internal boundary line segment, but less than the full length of the corresponding internal boundary line segment; (d) generating near internal boundary grid points between intersections of internal boundary intersected circles of adjacent geometry data points; (e) prioritizing the generated near internal boundary grid points by assigning weighed values with higher priority to the generated near internal boundary grid points which satisfy the near internal boundary grid spacing requirement; (f) performing unconstrained Delaunay triangulation on the prioritized near internal boundary grid points with the prioritized near internal boundary grid points being triangle vertices to form internal boundary grid line segments; (g) generating Voronoi cells to form unstructured gridding data with Voronoi cell edges corresponding to the internal boundary grid line segments; and (h) outputting the unconstrained gridding data. (a) refining the internal boundary geometry data to produce a partial set of geometry data points to represent the internal boundary geometry taking into account a near internal boundary grid spacing requirement, said refining comprises selecting a partial set of geometry data points where the distance between adjacent geometry data points satisfies the near internal boundary grid spacing requirements; (b) building lines by connecting adjacent geometry data points from the step above to form internal boundary line segments ; (c) constructing internal boundary intersected circles at each of the geometry data points, wherein the radius of each intersected circle is selected so that it is greater than a half-length of a corresponding internal boundary line segment, but less than the full length of the corresponding internal boundary line segment; (d) generating near internal boundary grid points between intersections of internal boundary intersected circles of adjacent geometry data points; (e) prioritizing the generated near internal boundary grid points by assigning weighed values with higher priority to the generated near internal boundary grid points which satisfy the near internal boundary grid spacing requirement; (f) performing unconstrained Delaunay triangulation on the prioritized near internal boundary grid points with the prioritized near internal boundary grid points being triangle vertices to form internal boundary grid line segments; (g) generating Voronoi cells to form unstructured gridding data with Voronoi cell edges corresponding to the internal boundary grid line segments; and (h) outputting the unconstrained gridding data. <EOS>