Derivative after numerical Integration

First I'll provide a symbolic workaround, and then I'll explain why your attempt doesn't work. The integral can be symbolically evaluated, like this:

ModelInternalEnergy[Td_, T_] := Evaluate[Simplify[ Alpha*T^4* Integrate[x^3/(Exp[x] - 1), {x, 0, Td/T}, Assumptions -> Td/T > 0] + ground]] 

which for reference gives

ground - 1/15 Alpha \[Pi]^4 T^4 + I Alpha \[Pi] T Td^3 - ( Alpha Td^4)/4 + Alpha T Td^3 Log[-1 + E^(Td/T)] + 3 Alpha T^2 Td^2 PolyLog[2, E^(Td/T)] - 6 Alpha T^3 Td PolyLog[3, E^(Td/T)] + 6 Alpha T^4 PolyLog[4, E^(Td/T)] 

You can then plot it, like this:

Plot[Re[ModelInternalEnergy[670, T] /. {Alpha -> 2.0, ground -> 3.0}], {T, 0, 10}, PlotPoints -> 30] 

You can then take the derivative symbolically, like this:

Simplify[Derivative[0, 1][ModelInternalEnergy][Td, T]] 

giving

-(4/15) Alpha \[Pi]^4 T^3 + I Alpha \[Pi] Td^3 + ( Alpha E^(Td/T) Td^4)/(T - E^(Td/T) T) + 3 Alpha Td^3 Log[1 - E^(Td/T)] + Alpha Td^3 Log[-1 + E^(Td/T)] + 12 Alpha T Td^2 PolyLog[2, E^(Td/T)] - 24 Alpha T^2 Td PolyLog[3, E^(Td/T)] + 24 Alpha T^3 PolyLog[4, E^(Td/T)] 

Your original attempt doesn't work because when you call Derivative[0, 1] on the NIntegrate expression, the NIntegrate expression is fed symbolic parameters T and Td, which makes it unhappy (NIntegrate needs to have purely numeric expressions in order to work).

In contrast, you can take derivatives of symbolic expressions, which is what the above does.

You may have also noticed that the plot used Re[...] of the model. This is because the symbolic expression has machine-epsilon-sized imaginary parts, which makes Plot unhappy, so Re is applied to remove those.

Also, note the Evaluate wrapper around the integral; this forces immediate evaluation of the expression, which is then SetDelayed to the result. := by default is a delayed operator, so without the Evaluate wrapper, it would evaluate the integral every time you call the function, which is not what we want.

If the model was too complicated for a symbolic solution, you could just manually define a numerical derivative in the T direction by finite differences (there are probably better options, though, so I'll leave that up to others to explain).

In V10,D can symbolically differentiate NIntegrate if it is Inactive. In the case that the integral can only be evaluated with NIntegrate the following defines an arbitrary Derivative of the OP's function numericalModelInternalEnergy.

ClearAll[numericalModelInternalEnergy, numericalModelHeatCapacity]; Block[{NIntegrate, Alpha, ground, T, Td}, integrand = Inactivate[Alpha*T^4*NIntegrate[x^3/(Exp[x] - 1), {x, 0, Td/T}] + ground]; numericalModelInternalEnergy[Td_?NumericQ, T_?NumericQ] = Activate[integrand]; Derivative[m_, n_][numericalModelInternalEnergy][Td0_?NumericQ, T0_?NumericQ] := Activate[D[integrand, {Td, m}, {T, n}] /. {Td -> Td0, T -> T0}]; ]; numericalModelHeatCapacity[Td_, T_] := Derivative[0, 1][numericalModelInternalEnergy][Td, T]; Block[{Alpha = 2.0, ground = 3.0}, Plot[{numericalModelInternalEnergy[670, T], 1000 numericalModelHeatCapacity[670, T]}, {T, 0, 1000}] ] 

The rule for the Derivative is stored as a sub-value:

SubValues[Derivative] (* {HoldPattern[Derivative[m_, n_][numericalModelInternalEnergy][ (Td0_)?NumericQ, (T0_)?NumericQ]] :> Activate[D[integrand, {Td, m}, {T, n}] /. {Td -> Td0, T -> T0}]} *) 

The definitions can be cleared with UnSet or by setting the SubValues directly:

Derivative[m_, n_][numericalModelInternalEnergy][Td0_?NumericQ, T0_?NumericQ] =.; 

or

SubValues[Derivative] = {} (* deletes all other user definitions, too *) 

Remark: The dependence of the definitions on the global symbols Alpha and ground is usually to be avoided. They could be made parameters of the functions, with ?NumericQ protection used for Td and T above.