In this work, the Precise Science Orbit (PSO) of the Gravity Field and Steady-State Ocean Circulation Explorer Mission (GOCE) satellite, acquired via the European Space Agency, served as the reference orbit for the testing various variants of the GOCE satellite orbit modeling. The GOCE satellite positions along the reduced-dynamic PSO orbit were treated as pseudo-observations in the satellite orbit improvement process in the least squares sense. This process was realized using dedicated extended Toruń Orbit Processor software which enabled determining corrections to the orbital initial conditions and the set of parameters, necessary to determine additional empirical accelerations of the satellite in the radial-track, along-track and cross-track directions. These piecewise accelerations were determined in successive equal intervals (pieces) into which the processed orbital arcs were divided. For modeling the accelerations, polynomials of different degrees were used. The obtained RMS differences between the improved orbits and the reference PSO orbit were determined for various orbital arc lengths up to 1-day. The best RMS of the fits for the 1-day arcs was in the range from 2.0 to 3.2 mm with significantly worst results for along-track direction. Based on the set of solutions determined, the number of orbital parameters for adopted accuracy thresholds and the upper limit of their number, that can be estimated depending on the length of the orbital arc, under given numerical conditions were obtained. The optimal form of the polynomial modeling of the estimated accelerations also depends partly on the length of the processed orbital arc. For shorter arcs (45 min and less), the second- third- and fourth-order polynomial gives the best results, while for longer arcs (90, 180, 360, 720 and 1440 min), zero- and first-degree polynomials are the most effective. A very promising solution with the RMS of the fit of 0.1 mm for a 1-min arc using the fourth-order polynomial was obtained in a perspective of the future use of the short arc approach to the fit of 1-day arcs. Additionally, solution variants with non-equal numbers of orbital pieces for the radial-along-and cross-track directions were obtained. In the case of the solution for an arc length of 1-day, the distribution of residuals and their statistics in the aforementioned radial-, along-and cross-track directions are presented.