By considering model comet nuclei with a wide range of sizes, prolate ellipsoidal shapes, spin axis orientations, and surface activity patterns, constraints have been placed on the nucleus properties of the primary Rosetta target, Comet 67P/Churyumov-Gerasimenko. This is done by requiring that the model bodies simultaneously reproduce the empirical nucleus rotational lightcurve, the water production rate as function of time, and non-gravitational changes (per apparition) of the orbital period (Delta P), longitude of perihelion (Delta omega tilde), and longitude of the ascending node (Delta Omega). Two different thermophysical models are used in order to calculate the water production rate and non-gravitational force vector due to nucleus outgassing of the model objects. By requiring that the nominal water production rate measurements are reproduced as well as possible, we find that the semi--major axis of the nucleus is close to 2.5 km, the nucleus axis ratio is approximately 1.4, while the spin axis argument is either 60+/-15 or 240+/-15 degrees. The spin axis obliquity can only be preliminary constrained, indicating retrograde rotation for the first argument value, and prograde rotation for the second suggested spin axis argument. A nucleus bulk density in the range 100-370 kg/m^3 is found for the nominal Delta P, while an upper limit of 500 kg/m^3 can be placed if the uncertainty in Delta P is considered. Both considered thermophysical models yield the same spin axis, size, shape, and density estimates. Alternatively, if calculated water production rates within an envelope around the measured data are considered, it is no longer possible to constrain the size, shape, and spin axis orientation of the nucleus, but an upper limit on the nucleus bulk density of 600 kg/m^3 is suggested.